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ABSTRACT 

We present a series of simulations of the fragmentation of a molecular cloud, leading to the 
formation of a cluster of protostellar cores. We use Gaussian initial conditions with a power 
spectrum P{k) oc fc -2 , assume an isothermal equation of state, and neglect turbulence and 
magnetic fields. The purpose of these simulations is to address a specific numerical problem 
O . called artificial fragmentation, that plagues simulations of cloud fragmentation. We argue that 

this is a serious problem that needs to be addressed, and that the only reasonable and practical 
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way to address it within the Smoothed Particle Hydrodynamics algorithm (SPH) is to use a 
CI relatively new technique called particle splitting. 

We performed three simulations, in which we allow iV gcn = 0, 1, and 2 levels of particle split- 
ting. All simulations start up with 64 3 SPH particles, but their effective resolutions correspond 
to 64 3 , 128 3 , and 256 3 particles, respectively. The third simulation properly resolves the Jeans 
mass throughout the entire system, at all times, thus preventing artificial fragmentation. 
■ The high resolution of our simulations results in the formation of a large number of protostellar 

cores, nearly 3000 for the largest simulation. This greatly exceeds the typical number of cores 
(~ 60) formed in previous simulations, and enabled us to discover various processes that affect 
the growth of the cores and the evolution of the cluster. 

The evolution of the cloud follows four distinct phases, or regimes. Initially, during the growth 
regime, the cloud evolves into a network of intersecting filaments. After roughly one dynamical 
time, core formation starts inside the dense gaseous fragments located at the intersection of the 
filaments, and to a lesser extent inside the filaments themselves, and the cloud enters the collapse 
regime. During this regime about 50% of the gas, essentially the gas that started up in overdense 
regions, is converted into cores. Competitive accretion is the main process that controls the 
mass evolution of cores, but we discovered that this process occurs locally, within each dense 
fragment. During the following accretion regime, most of the remaining gas does not form new 
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cores, but rather accretes onto the existing cores. Eventually each gaseous fragment has turned 
into a subcluster of cores, and these subclusters later merge to form the final cluster. The gas 
left in the system has become negligible, and the system has reached the N-body regime, in which 
the dynamics of the cluster is governed by N-body dynamics. 

The final mass distribution of cores has a lognormal distribution, whose mean value is 
resolution-dependent; the distribution shifts down in mass as the resolution improves. The width 
of the distribution is about 1.5 (e.g., a factor of 30 in the mass), and the low-mass edge of the 
distribution corresponds to the lowest core mass that the code can resolve. This result differs 
from previous claims of a relationship between the mean of the distribution and the initial Jeans 
mass. 

Subject headings: hydrodynamics — ISM: clouds — methods: numerical — stars: formation 
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1. Introduction 

1.1. Clustered Star Formation 

An understanding of star formation is pivotal to understanding both the origin of galaxies and the 
origin of planetary systems. A crucial step in the formation of galaxies is the formation of stars. Intensive 
archaeology of chemical abundances has uncovered the outline of the star formation history of our own galaxy 
(e.g., McWilliam 1997), and look-back studies are beginning to provide more direct information about star 
formation beyond z ~ 1. It now appears that the star formation rate was flat to z > 3 (Madau 1999; Smail 
et al. 2000) or even increasing up the highest observed redshifts (Lanzetta et al. 2002). At least half the 
extragalactic background light lies in the far-infrared to submillimeter region (Puget et al. 1996; Hauser 
et al. 1998), suggesting that half the star formation in the history of the Universe has occurred in dusty 
regions. Submillimeter galaxies, at a mean z = 2.3 (Chapman et al. 2005), are forming stars at prodigious 
rates (up to 1000 solar masses/year), in dusty, molecule-rich environments (e.g., Barger et al. 1998; Cimatti 
et al. 1998; Cox et al. 2002). Such objects have now been seen even at z — 6.4 (Walter et al. 2004). To 
understand the nature of this process, we must develop a better understanding of massive star formation in 
similar environments more amenable to detailed study, the dense regions of molecular clouds in our Galaxy. 

In a similar way, the study of the origin of planetary systems rests on a deeper understanding of the 
role of disks in star formation. It is now clear that disks surround most young stars (Beckwith et al. 1990), 
and we are beginning to study the properties and evolution of these disks (e.g., Mundy, Looney, & Welch 
2000) . These studies are fundamental to understanding the frequency and variety of planetary systems that 
current searches are revealing (e.g., Cumming, Marcy, & Butler 1999; Jorissen, Mayor, & Udry 2001; Marcy 
et al. 2005). 

Over the last two decades, we have made substantial progress in understanding how low-mass stars 
form in relative isolation. Spurred by the elaboration of an elegant theoretical paradigm, beginning with Shu 
(1977) and culminating in the influential review by Shu, Adams, & Lizano (1987), observers have developed 
the capability to test predictions of theoretical models. We have a "standard model" and a number of 
variations, each of which makes predictions that can be tested by observations. In particular, the form of the 
density and velocity field in the envelope around the forming protostar differs in the different models, and 
observers are beginning to be able to distinguish these observationally (see reviews by Evans 1999; Myers, 
Evans, & Ohashi 2000; Andre, Ward-Thompson, & Barsony 2000). Studies of dust continuum emission with 
instruments like SCUBA is providing a valuable new probe of the density distribution, while molecular line 
profiles probe the kinematics (e.g., Zhou et al. 1993; Zhou & Evans 1994; Myers, Evans, & Ohashi 2000). In 
addition, predictions that disks would form on scales of 1 — 100 AU around forming stars (e.g., Terebey, Shu, 
& Cassen 1984) have been verified (Beckwith et al. 1990). In the study of low-mass, isolated star formation, 
theory has clearly revealed the path to the observers. 

However, most stars probably form not in isolation, but in clustered environments (Elmegreen 1985; 
Clarke, Bonnell, & Hillenbrand 1999; Evans 1999; Pudritz & Fiege 1999; Elmegreen et al. 2000; Lada & 
Lada 2003). Massive stars seem to form exclusively in these environments, but the full spectrum of stars and 
sub-stellar objects forms in clusters (Lada & Lada 1995). In addition, the bursts of star formation seen in 
distant galaxies are clearly related to the formation of massive stars in clusters. Consequently, understanding 
them is crucial to understanding galaxy formation. Disks around forming stars have been seen in clustered 
environments (e.g., O'Dell & Wen 1994; Strom 1995), but interactions among proximate star/disk systems 
may affect the mass distribution of these disks (Eisner & Carpenter 2003). 
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In the area of massive, clustered star formation, extensive observations exist (e.g., Churchwcll 1993; 
Kurtz et al. 2000), but theory is very underdeveloped. We know that the typical conditions in dense regions 
of high mass are quite different from those in dense regions of low mass (e.g., Plume et al. 1997; Shirley et 
al. 2003b). The temperatures and densities are higher, and the masses of gas are sufficient to form many 
stars (up to 10 4 M Q ). The linewidths of molecular lines are much greater and deviate from the linewidth-size 
relations known from lower-mass regions and regions not forming stars (Shirley et al. 2003a). The overall 
density distribution in these high-mass regions seems well approximated by a power law n{r) = rif(r/rf)~ p 7 
with exponent similar to that of low-mass regions (p ~ 1.8) (Shirley et al. 2003b) but densities (n/) 100 
times higher (Evans et al. 2002; Mueller et al. 2002). These wide lines are highly supersonic, and turbulence 
must play a major role. Some dense regions seem to be more fragmented (Wang et al. 1993; Lada, Evans, & 
Falgarone 1997), as might be expected if clusters are forming. 

In these conditions, direct extension of theories of low mass star formation may run into problems. 
Recent work on the formation of massive stars delineates the theoretical framework (McKee & Tan 2002, 
2003) of formation in turbulent dense regions. This work predicts the overall properties of the regions, but the 
details of the formation and its observational manifestations remain to be understood. The close proximity of 
other clumps within the overall dense region will perturb the density and velocity fields around a given clump, 
making the kinds of tests applied to isolated regions less meaningful (and very difficult observationally) . 
Theories that make detailed predictions of observables are needed. Statistical measures must be compared 
between theory and observation, rather than any specific realization of a simulation because stochasticity 
is inevitable. For example, the distribution of clump masses and velocities as a function of time could be 
predicted and compared to observations. If one can follow the process with sufficient dynamic range, the 
distribution of clump angular momenta could be used to predict things like the frequency of binaries and 
disks, if supplemented by more detailed calculations of the subsequent evolution of individual structures once 
they get small enough that their internal dynamical time becomes less than the interaction time. Indeed, 
predictions of when that point occurs are needed. Larger scale correlations can be studied, such as the 
tendency toward alignment of angular momentum vectors or magnetic fields. 

A number of simulations of the fragmentation of a molecular cloud to form a cluster have been performed 
(e.g., Larson 1978; Keto, Lattanzio, & Monaghan 1991; Bonnell et al. 1997; Bromm, Coppi, & Larson 1999; 
Padoan & Nordlund 2002; Tilley & Pudritz 2004; see the recent review articles by Larson 2003; Mac Low & 
Klessen 2004). Two research groups have been particularly prominent in recent years. Klessen, Burkert, and 
their collaborators (Klessen, Burkert, & Bate 1998; Klessen & Burkert 2000, 2001; Klessen, Heitsch, & Mac 
Low 2000; Klessen 2001a, b; Schmeja & Klessen 2004; Jappsen & Klessen 2004, hereafter collectively KB; 
Li, Klessen, & Mac Low 2003; Jappsen et al. 2005), have used a Smoothed Particle Hydrodynamics (SPH) 
code to follow the evolution of a region with many times the Jeans mass. They were able to reproduce 
the observed distribution of clump masses, but found that bound units developed a steeper distribution, 
more similar to that of stars. This result meshed nicely with recent observations showing that the densest 
structures in several clouds also have a steeper mass function for masses above 0.5M Q (Motte, Andre, & 
Neri 1998; Testi & Sargent 1998). Bonnell, Bate, and their collaborators constitute the second group (Bate, 
Bonnell, & Bromm 2002a,b, 2003; Bate & Bonnell 2005; Dale et al. 2005; Dobbs, Bonnell, & Clark 2005). 
They also used SPH, but have focused mostly on smaller, denser systems for which the gas is optically thick 
at high densities and the assumption of isothermality breaks down. 

An important issue in such simulations is the phenomenon called artificial fragmentation, which can 
lead to severe numerical problems. In the next subsection, we describe the approach used by the various 
groups for dealing with this problem. 
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1.2. Artificial Fragmentation and the Jeans Criterion 

Truelove ct al. (1997) and Boss (1998) have shown that a minimum requirement of any grid-based 
simulation of a fragmenting cloud is that the algorithm can resolve the Jeans mass Mj (we shall refer to 
this as the Jeans criterion). The maximum mass inside a cell must be smaller than ~ 1/64 of the Jeans 
mass in order to prevent a spurious, resolution-dependent effect they called artificial fragmentation. This 
effect, if present, can invalidate the results of star-formation simulations, by producing initial mass functions 
and accretion histories that are totally wrong. Bate & Burkcrt (1997) derived a different, but equivalent 
criterion for SPH. The total mass contained inside the zone of influence of a particle must be less than 
about twice the Jeans mass in order to prevent artificial fragmentation. As Klein ct al. (1999) pointed out, 
this requirement poses a serious problem for isothermal simulations performed with SPH. The Jeans mass 
varies as Mj oc T 3 ! 2 p~ x l 2 . Hence, in isothermal clouds, Mj oc p^ 1 ^ 2 , so the Jeans mass decreases as collapse 
proceeds. Since the smoothing lengths h are adjusted in such a way that the mass inside the zone of influence 
of every particle remains constant, the Jeans mass will eventually be underresolved as the density increases. 

There is, however, a physical lower limit to the mass of fragments, simply because the isothermal 
approximation breaks down at sufficiently high densities. Burkert, Bate, & Bodenheimer (1997), Bate 
(1998), and Klein ct al. (1999) extend their simulations into the high density regime by using a barotropic 
equation of state, which is isothermal below a certain critical density p C 2 , and adiabatic above it (the actual 
equation of state is significantly more complex [see, e.g., Scalo et al. 1998], but the barotropic form is a 
convenient approximation). In the adiabatic regime, P oc p 5 ^ 3 , hence T oc p 2 ! 3 and therefore Mj oc p 1 ! 2 . 
There is therefore a minimum Jeans mass, corresponding to the critical density p C 2. As long as the Jeans 
criterion is satisfied at that density, it will be satisfied in the entire system, at all times. 

This still poses a problem. The isothermal approximation is valid for densities up to 10 10 cm -3 (Klessen 
& Burkert 2000). For a region with mean density 10 2 cm~ 3 , the range of 10 8 in density corresponds to a 
range of 10 4 in Jeans mass. Since the smallest Jeans mass must contain ~ 100 particles to satisfy the Jeans 
criterion, the total number of particles in the simulation must be at least 1 million, and this would satisfy 
the Jeans criterion only marginally. 

Klessen, Burkert, and their collaborators do not address the problem of artificial fragmentation. They 
are aware of the problem, but believe that this problem does not affect their conclusions significantly. Bate, 
Bonnell, and their collaborators are very concerned with artificial fragmentation. They avoid the problem by 
simulating low-mass systems (M ~ 50M Q ), so that their resolution is sufficiently large to reach the regime 
where the gas becomes adiabatic. They can then use a barotropic equation of state. 

In this paper, we consider an alternative approach that can simulate high-mass systems, in the isothermal 
regime, while still solving the artificial fragmentation problem. This is achieved by using particle splitting, 
a relatively new and very promising technique (Kitsionas & Whitworth 2002). The basic idea consists of 
starting the simulations with a manageable number of SPH particles, and then refining the mass resolution 
locally in regions where additional resolution is needed to satisfy the Jeans criterion. Original particles 
are automatically replaced (or "split" ) by a more finely-spaced set of smaller-mass particles wherever extra 
resolution is required. This can be seen as the Lagrangian counterpart of the Adaptive Mesh Refinement 
techniques used in Eulerian, grid-based algorithm. In our implementation of particle splitting, SPH particles 
split into 8 equal-mass particles when the Jeans criterion is locally violated. For an isothermal gas, this 
results in a new generation of particle splitting every time the density increases by an additional factor of 
64. 
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1.3. Objectives 

The primary goal of this ongoing project is to study the effect of feedback from clustered star formation 
on the evolution of the ISM. However, the issue of feedback will not be addressed in this first paper, for the 
following reason: The existence of artificial fragmentation casts a huge shadow over all SPH simulations of 
cloud fragmentation that assume an isothermal equation of state, raising doubts about the validity of such 
simulations. We feel that this problem is far too important to be ignored, and must be addressed first, and 
successfully, before we even consider implementing additional physical processes into the algorithm. Solving 
the problem of artificial fragmentation is an essential first step. In this paper, we address this problem using 
a SPH algorithm which combines self-gravity, hydrodynamics, particle splitting, and sink particles. This 
algorithm is described in §2 below. 

The main objectives of this paper are the following: 

1. Test the feasibility of particle splitting, and investigate, both analytically and numerically, the interplay 
between particle splitting, sink particles, and the Jeans criterion. Using particle splitting enables us to start 
up simulations with a small number of particles for a given resolution. However, the feasibility of this 
approach depends critically on the efficiency of particle removal by sink formation and accretion onto sinks. 
If a large number of particles split before the formation and growth of sinks becomes important, the total 
number of particles in the simulation might become too large to be manageable. It is not obvious a priori 
that sink formation and growth will remove particles sufficiently rapidly to offset the increase in particle 
number resulting from splitting. One of the main goals of this paper is to determine if the peak number of 
particles in such simulations remains manageable. 

2. Perform a convergence study. We can ascribe to each simulation with particle splitting an "effective 
particle number," which is the number of particles a simulation without particle splitting would need to 
achieve the same resolution in dense regions. The three simulations presented in this paper start with 64 3 
particles, and allow for N gen = 0, 1, and 2 generations of particle splitting, respectively. The splitting factor 
is /split = 8, meaning that when a particle splits, it is replaced by 8 particles, each having 1/8 of the mass of 
the parent particle. Hence, the effective particle numbers of the three simulations are 64 3 , 128 3 , and 256 3 , 
respectively. Since we are using identical initial conditions, these three simulations taken together constitute 
a convergence study, the largest one ever performed for such simulations. 

3. Perform the largest simulation of this kind ever done, in terms of effective number of SPH particles 
or number of proto stellar cores formed. For our largest simulation (A gcn = 2), the effective particle number 
is 256 3 = 16,777,216, about 33 times the largest number of particles used by KB (500,000). Our smallest 
simulation (N gcn = 0, that is, no particle splitting) uses 64 3 = 262, 144 particles, which is comparable to the 
largest isothermal simulations of KB. Furthermore, we start the simulations with Nj = 500 Jeans masses 
instead of KB's 222 Jeans masses. As a result, we will form a much larger number of protostellar cores. This 
has three advantages: First, our determination of the initial mass function of protostellar cores will be more 
accurate. Second, it will enable us to study the mass assembly history and final structure of the cluster in 
more detail. And third, forming a larger number of cores might lead to the discovery of some interesting 
processes in the evolution of the cluster, that would not occur in a cluster with much fewer cores; as we shall 
see, this is indeed the case. 
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2. THE NUMERICAL ALGORITHM 



2.1. Basic Equations 



The evolution of a self-gravitating gas is described by the conservation equations for mass, momentum, 
and energy, coupled with the Poisson equation and the equation of state, 



^+V-(pv) = 0, (1) 
<9v VP 

7*+(vV)v = -— -V*. (2) 

g +v .Ve=-PV.v + H, (3) 

V 2 = 4irG(p - p) , (4) 

P = f(p,e), (5) 



where p is the density, P is the pressure, e is the specific internal energy, v is the velocity, (j> is the gravitational 
potential, p is the mean density, and T and A are the radiative heating and cooling rates, respectively. 
Equation (4) requires some explanation. To prevent the overall collapse of the cloud, we assume that the 
cloud is essentially infinite, and use periodic boundary conditions. However, a periodic gravitational potential 
(f> is only possible if the total mass of the system vanishes. By adding the term — p in equation (4), the effective 
mass, defined as the integral of the source term p — p over the computational volume, does vanish, and a 
periodic solution for (f> becomes possible. A valid interpretation of equation (4) is that the term — p accounts 
for whichever process makes the cloud globally stable, while it is the fluctuation p — p that make the cloud 
locally unstable. Klessen, Heitsch, & Mac Low (2000) and Mac Low & Klessen (2004) have suggested that 
supersonic turbulence might explain the global stability of clouds. 

It is important to realize that equation (4) is still physically correct. Since the cloud is assumed to be 
infinite, the term — p represents a uniform, negative density background extending to infinity in all directions, 
and such component cannot, by symmetry, exert any force in any direction on a mass element. For a more 
formal description, we refer the reader to Alccian & Leorat (1998), and references therein. 

In this paper, the set of equations we are using is significantly simpler. First, since we use Smoothed 
Particle Hydrodynamics (SPH), a Lagrangian, particle-based algorithm, mass is automatically conserved, 
and we can ignore the continuity equation (1). Second, we assume that the gas is isothermal. The specific 
internal energy e is therefore constant in space and time, and equation (3) can be ignored as well. The 
equation of state, equation (5), becomes 



where c s (e) is the sound speed, which is constant in an isothermal gas, and 7 is the polytropic constant 1 . 



1 Thc relationships between the concepts of polytropic constant, polytropic equation of state, and isothcrmality are often 
reported incorrectly. The polytropic constant 7 is the ratio of the specific heat capacities at constant pressure and constant 
volume. A polytropic equation of state has the form P oc p 7 , but such an equation is valid only if the entropy of the gas is 
constant both in space and time. It is often said that 7 = 1 for an isothermal gas, but this is incorrect in general. In our 
simulations, 7 = 5/3, the equation of state is not polytropic, and it is the coupling between the gas and a background radiation 
field that makes the gas isothermal. 
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The system of equations (l)-(5) reduces to 



_^ V p-V0, (7) 
at 7p 

V 2 = A7rG(p - p) , (8) 



where d/dt = d/dt + v • V is the Lagrangian time derivative. 



2.2. The SPH/P 3 M Algorithm 

We use a hybrid gravity/hydrodynamics algorithm. The gravitational forces are computed using a 
P 3 M algorithm (Hockney & Eastwood 1981), while the gasdynamical equations are solved using the SPH 
algorithm (see Monaghan 1992, and references therein). This hybrid algorithm was originally introduced by 
Evrard (1988), though we have developed our own version. Our code is actually an Adaptive SPH (ASPH) 
code (Shapiro et al. 1996; Owen et al. 1998). The ASPH algorithm uses anisotropic smoothing kernels, and 
a special treatment of artificial viscosity that prevents spurious preheating of low-density regions. However, 
in this paper, we consider an isothermal equation of state, so preheating is not an issue. We decided to use 
isotropic smoothing kernels for now (that is, using the Adaptive SPH code as a standard SPH code). We 
will consider anisotropic smoothing kernels in future work. 

We have modified our original algorithm to include a treatment of particle splitting and sink particles. 
The implementation of these features is described below. 



2.3. Particle Splitting 

We have modified our original SPH algorithm to include particle splitting and sink particles. Particle 
splitting is implemented as follows. We assume that there is a minimum number of particles nj >a .f. that a 
Jeans mass must contain in order to be properly resolved and not undergo artificial fragmentation. This 
gives us the Jeans criterion. Each particle i is required to satisfy the condition 

> nj, a .f. , (9) 

mi 



where Mj is the Jeans mass, which is a function of the density p and specific internal energy e, Mj(pi 1 ei) 
is the value of Mj evaluated at the location of particle i (see §4.3 below), and mi is the mass of particle i. 
When this condition is violated, the algorithm responds by splitting particle i into / sp iit particles of mass 
m = nii/ f 'spin- These new particles might themselves be split later as Mj keeps decreasing. For instance, in 
the case of an isothermal calculation (e = const), Mj oc p~ 1/>2 , and there is a series of particular densities 
p = ap, a/0.f s 2 plit , apf^ plit , . . ., (a = const) at which particle splitting will occur. 

When a particle is split, the algorithm must determine the positions, velocities, masses, and smoothing 
lengths of the daughter particles 2 . In the original approach of Kitsionas & Whitworth (2002), particles are 
split into a sphere of 13 daughter particles forming a compact lattice. We use a different approach, in which 



2 and also the specific internal energies, if the gas is not isothermal. 
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particles are split into / sp iit = 8 particles located on the vertices of a cube. This approach will be easier to 
generalize to Adaptive SPH later. We create the daughter particles as follows: Consider a particle i, with 
position r i: velocity Vj, mass m„ and smoothing length hi, that violates the condition (9), and therefore 
needs to be split. If (Ar)i is the mean particle spacing in the vicinity of particle i, the 8 daughter particles 
will be located at 

, (Mi 



±i 
±i 



(10) 



so that the spacing between daughter particles will be (Ar)j/2. 3 In the initial conditions, the smoothing 
lengths hi of the particles are initialized to be a multiple of the mean particle spacing: ft-unit = ^Ar, where 
£ 2 is a constant (using the notation of Shapiro et al. 1996). Then, as the calculation proceeds, the algorithm 
evolves the smoothing lengths hi is such a way that this relation is (roughly) maintained locally. Hence, the 
local particle spacing can be estimated using (Ar)j = hi/£ 2 - Equation (10) reduces to 



4/ 2 



"±1" 
±1 
±1 



(11) 



In our simulations, we generated initial conditions using £ 2 = 2, to provide a sufficient number of neighbors. 

We set the velocity of the daughter particles equal to Vj. We could use a more sophisticated approach 
that would take the local velocity gradient into account, but this is not really necessary. Particles tend to 
readjust themselves in one time step, erasing any velocity fluctuation at scales smaller that the smoothing 
length. Finally, we set the masses of the daughter particles equal to rrii/8, and their smoothing lengths equal 
to hi/2. 



2.4. Sink Particles 

Sink particles must be used in simulations of cloud fragmentation to reduce the timesteps, if we hope 
to determine the initial mass fraction of collapsed fragments. It is not sufficient to use an algorithm with 
individual timesteps, because as the gas fragments and collapses, most SPH gas particles end up in dense 
regions with short dynamical times. By replacing each dense clump of gas particles by a single object, called 
a sink particle (Bate, Bonnell, & Price 1995; Bromm, Coppi, & Larson 2002), we can eliminate this problem 
and increase the speed of the algorithm tremendously. In our implementation of sink particles, we use the 
method of Bate, Bonnell, & Price (1995), and we refer the reader to that paper for details. 



2.4-1- Creation of Sinks 

Our algorithm uses individual timesteps. Each particle i is given a timestep AU — (Ai)b a sic/2™, where 
(At)basic is the "basic timestep," and n identifies the timestep bin where the particle belongs. Sink particles 
are created only at the end of the basic timesteps, when all the particles are in sync. This greatly simplifies 



3 Note: because of the periodic boundary conditions, some daughter particles are wrapped around the computational box if 
particle i is located less than (Ar); /4 away from the edge of the box. 
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the implementation. Gas particles are replaced by sinks when they reach a threshold density p c . This density 
threshold is a numerical parameter that does not have any physical meaning, except for the fact that the 
assumption of isothermality must be physically valid at all densities below p c . The smallest collapsed objects 
that can form in the simulations will have a mass equal to (Mj) c the Jeans mass at the density p c . Hence, 
p c fixes the mass resolution of the algorithm, much in the same way as the softening length fixes its length 
resolution. 

To create sinks, the algorithm identifies all particles whose density 4 exceeds the threshold density p c 
for sink creation. These particles are then sorted in decreasing order of the density. Starting with the 
densest particle, the algorithm finds all n acc particles within a fixed radius r acc of that densest particle, 
where r acc is called the accretion radius. An important issue is how to determine the appropriate value for 
r acc , or equivalcntly, how many particles n acc should be turned into a sink. We shall assume that it takes a 
minimum number of particles, nj jm i n , to properly resolve a Jeans mass, and we will set that number equal 
to the number of particles nj >a .f. that a Jeans mass must contain to prevent artificial fragmentation. This 
makes sense, since the negative consequence of underresolving the Jeans mass is precisely to cause artificial 
fragmentation. The number of particles n acc must also exceed the number of particles (nj) c inside a Jeans 
mass at the threshold density p c , otherwise sub- Jeans mass objects would be turned into sinks. Hence, the 
condition is 

n acc > max [nj, min , (nj) c ] . (12) 

This ensures that every fragment replaced by a sink particle (i) is properly resolved, and (ii) contains at 
least a Jeans mass. We adjust r acc such that n acc « max[nj jm i n , {nj) c \ at the threshold density p c . If the 
values of n acc are systematically smaller than max [nj jm i n , [nj) c ], this would imply that we chose a value 
of r acc that was too small, and sinks would then form at densities much larger than p c . Conversely, if the 
values of n acc are systematically larger than max [nj jm i n , [nj) c ], this would indicate a value of r acc that was 
too large, effectively reducing the mass resolution of the algorithm, by forming objects that are too massive. 
The actual determination of r acc is described below. 

To create a sink particle, a second condition must be satisfied, 

"lib* 1 - (13 » 



where E t h and E gI are the thermal and gravitational energies of the particles inside the accretion radius, 
respectively. Bate, Bonnell, & Price (1995) use a < 0.5 as a criterion. However, using equation (13) is more 
consistent with the fact that we are trying to turn Jeans-mass clumps into sinks, since by definition a = 1 
for a uniform sphere of mass M = Mj. If both conditions are satisfied, the particles inside r acc are removed 
and replaced by a sink particle, which inherits the properties of the parent particles (center of mass position 
and velocity, total mass, and total angular momentum). The algorithm then selects the next densest particle 
still available (some particles with p > p c might have been incorporated into sinks created around denser 
particles), and the process is repeated until all particles with p > p c have been considered. 



4 We use the expression "density of a particle" for convenience; strictly speaking, SPH particles do not carry a density, hence 
the correct expression is "the density of the gas at the location of a particle." 



-10- 



2.4-2. Mass Evolution of Sinks 

Two physical processes can increase the mass of protostellar cores: mergers with other cores, and 
accretion of gas by cores. In the algorithm, these processes correspond to merging of sinks, and accretion of 
SPH gas particles by sinks, respectively. 

The merging of sinks is an interesting issue. Allowing sinks to merge would clearly affect the final mass 
distribution of protostellar cores. Unfortunately, no known prescription exists for implementing sink merging 
into a SPH algorithm. We intend to investigate this issue in a separate paper, but for now we will make the 
usual assumption that sink particles do not merge, as in Bate, Bonnell, & Price (1995), KB, and others. 

Accretion of gas particles onto sinks is also performed at the end of each basic timestep. The algorithm 
checks for gas particles that are located inside the accretion radius r acc of a sink particle. These particles are 
then accreted by the sinks, and the sink properties (position, velocity, mass, angular momentum) are updated 
accordingly. If a gas particle is within the accretion radius of several sinks, we compute the total energy of 
the (particle + sink) systems, and choose the sink for which this energy is the smallest. Unlike Bate, Bonnell, 
& Price (1995), we do not require that the energy be negative. Such a criterion would fail in general, because 
in these simulations the gas tends to fall onto sinks at supersonic speeds that exceed the escape velocity from 
the sink. This happens when gas particles are accelerated by large mass concentrations, that might contain 
several sinks. A gas particle approaching a cluster of sinks will be accelerated by the whole cluster, and will 
acquire a velocity comparable to the escape velocity from the whole cluster. But this velocity will always 
exceed the escape velocity from the particular sink that will accrete that particle. Hence, the gas particle 
will not be gravitationally bound to the sink onto which it accretes. This is not a concern, because the 
physical process responsible for the final stage of accretion is not gravitational capture, but rather physical 
collision between a gas particle and a sink. In a simulation without sinks, a gas particle approaching a dense 
clump at supersonic speed would be decelerated down to subsonic speeds by the artificial viscosity, resulting 
in the conversion of kinetic energy into thermal energy, that would then be radiated away since the gas is 
assumed to be isothermal. But this process would occur at scales smaller than r acc , and therefore cannot 
occur when clumps are being replaced by sinks. By allowing large-velocity particles to accrete onto sinks, 
we are essentially putting-in the subgrid physics of collision, viscous deceleration, and radiative dissipation 
by clumps smaller than r acc . 

When using sinks, boundary conditions must be applied at the accretion radius r acc , otherwise the SPH 
calculation of the density, pressure forces, and viscous forces on gas particles located immediately outside 
face would be incorrect, resulting in spurious effects. This is discussed in great detail in Bate, Bonnell, & 
Price (1995), and also in Bromm, Coppi, & Larson (2002). Our implementation of boundary conditions 
follows the description of Bate, Bonnell, & Price (1995). Interestingly, in our simulations, the effect of 
boundary conditions turns out to be quite small. We found that most gas particles fall radially toward sinks, 
and collide with them at supersonic speeds. As Bate, Bonnell, & Price (1995) point out, viscous boundary 
conditions have a negligible effect for particles falling radially, while pressure boundary conditions have a 
negligible effect for particles moving supersonically. Our findings are consistent with that claim. 

3. CLOUD FRAGMENTATION 

In this section, we discuss some aspects of the numerical simulations of cloud fragmentation, including 
the interplay between resolution, artificial fragmentation, particle splitting, and sink pyarticles. For the sake 
of the discussion, we consider a particular simulation of KB, which we call the "benchmark simulation." For 
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this particular simulation, the number of particles is N — 200, 000, the initial number of Jeans masses is 
Nj = 222, the mean density, in computational units, is p — 1/8, and the threshold density for creating sinks 
is p c = 40, 000p. For the number of particles necessary to prevent artificial fragmentation, we assume nj >a .f. « 
100 (this is a guess; no value of nj^.f. is provided by KB). Notice that this is one of the highest-resolution 
simulation performed by KB. They have performed a higher-resolution simulation with N = 500, 000, and 
by lowering the threshold density and the initial number of Jeans masses, they were able to increase their 
mass resolution significantly. Because of the reduced number of Jeans masses in this simulation, the more 
suitable calculation for us to use as a benchmark is the one with N — 200, 000. 



We can estimate the effect of particle splitting on the mass resolution of the algorithm using a simple 
calculation. We assume that the system has a total mass M to t, and the simulation starts up with N equal- 
mass particles of mass m = M tot /N. Because the initial density fluctuations are small (p; n it ~ p), and the 
gas is isothermal, the initial Jeans mass Mj^it is essentially constant in space. The initial number of Jeans 
masses in the system is then given by Nj = M tot /Mj^ n i t , and the initial number of particles in each Jeans 
mass is given by n Jiin it = N/Nj = M JMi t/m. 

In the initial state, the system must satisfy the following conditions: 



The first condition states that the Jeans mass is properly resolved initially, thus preventing immediate 
artificial fragmentation, while the second one states that the initial density fluctuations are too small to 
immediately trigger the formation of sink particles. As the system evolves, the density p increases in some 
regions, and since nj oc Mj oc p~ 1//2 for an isothermal gas, nj decreases in the same regions. Eventually 
the conditions (14) and (15) will both be violated in dense regions. If condition (15) is violated first, a 
sink particle will be created, and this will prevent artificial fragmentation from happening. If instead the 
condition (14) is violated first, then artificial fragmentation will occur, but as long as the fragments remain 
close to one another, they might eventually get lumped together into a single sink particle once that sink is 
created. The creation of a sink particle requires an increase in density by a factor / s i n k = p c /p~, while violating 
the Jeans criterion requires an increase in density by a factor / a .f. = (nj.init/nj.a.f.) 2 = {N /Njnj^.i.) . For 
the benchmark simulation, we get / s i n k = 40, 000 and / a .f. — 81. Hence, the Jeans criterion will be violated 
after an increase in density by a factor of 81, long before sink particles are created. Artificial fragmentation 
must be prevalent in the isothermal simulations of KB, but these fragments might get re-aggregated when 
sink particles are created. 5 . 

There are four possible solutions to the problem of artificial fragmentation. The most obvious solution 
(1) consists of reducing p c down to 81p or less, so that fragments would turn into sink particles before they 
are dense enough to experience artificial fragmentation. The obvious drawback of this approach is that 
the simulation would completely miss the late-stage evolution of dense fragments. If we keep the threshold 
density at a value p c = 40, 000p, we must then find a way to increase the number of particles per Jeans 
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5 Actually, the limited resolution of the gravity solver might help to prevent artificial fragmentation in the simulations of KB 
(Klessen 2002). 



mass. We could (2) increase the initial Jeans mass Mj^it by a factor of (40, 000/81) 1 / 2 ~ 22, so that by 
the time the density reaches p c , the Jeans mass would still contain enough particles to be properly resolved. 
However, for a fixed initial number TV of particles, this would reduce the initial number of Jeans masses in 
the system from Nj = 222 down to Nj — 10, leading to poor statistics, as very few cores would form. 

If we keep both p c and Nj fixed, we must then reduce the particle mass by a factor of at least 
(40,000/Sl) 1 / 2 to insure sufficient resolution. If (3) this is done over the entire computational volume, 
the number of particles would then increase from 200,000 to 4,400,000. This is clearly a brute-force ap- 
proach, that would result in a tremendous increase in computational time. The better approach (4) consists 
of reducing the particle mass locally, in dense regions where the Jeans mass is small. This solution can be 
achieved dynamically, using particle splitting. 

After splitting, the mass per particle is given by 

m = m init /~^ on , (16) 



where mj n i t = M to t/N is the initial mass of the particles, / sp iit is the splitting factor, and A gon is the 
maximum number of "generations," that is, the maximum number of splittings a particle can experience. 
At the threshold density p = p c , the Jeans mass is given by 



We eliminate mi n i t in equation (17) using equation (16). The Jeans criterion, (Mj) c /m > nj ia .f. becomes 



Nj\p c 



1/2 

/ B pST > nj,a.f. • (18) 



This relates the initial number of Jeans masses Nj to the number of splitting generations N gcn necessary to 
prevent artificial fragmentation. This expression is useful for generating initial conditions (see §4.2 below). 



3.2. Scenarios 



The evolution and fragmentation of the cloud can follow three different paths, which we call scenarios. 
The particular scenario followed is determined by two parameters: the maximum number of splitting gener- 
ations A gen and the accretion radius r acc , which determines, among other things, the number of gas particles 
that are lumped together when a sink is created. Figure 1 illustrates the three different scenarios. Consider 
first a simulation with an insufficient number of splitting generations, or no splitting at all (A gcn too small). 
As the system evolves, the density increases and the Jeans mass decreases in overdense regions. Eventually, 
the cloud fragments into Jeans-mass clumps, which themselves fragment as the density increases and Mj 
decreases, until the density reaches p a .f. = f a .f.p. At that point, artificial fragmentation occurs, leading to 
fragments of mass M < Mj. As p increases and Mj decreases, the mass of these "artificial fragments" might 
eventually end up above Mj, but further artificial fragmentation will occur, bringing the mass per fragment 
below Mj again. Finally, when the density reaches the threshold value p c , sink particles are created. KB 
suggest setting the accretion radius r acc equal to the Jeans length (Rj) c at the density p — p c . Each sink 
particle will then contain about one Jeans mass, implying that several sub- Jeans-mass fragments will be 
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Fig. 1. — Schematic illustration of the various scenarios. Open circles represent fragments composed of SPH 
particles. Solid dots represent sink particles. From left to right, time increases, local density increases, and 
Jeans mass decreases. 



lumped together into the same sink particle, possibly nullifying the effects of artificial fragmentation. This 
is Scenario I. 

One concern with this scenario is that the Jeans mass (Mj) c at the threshold density p c might be greatly 
underresolved (whether or not artificial fragmentation occurred) , and replacing it by a sink particle might be 
inappropriate. For instance, the benchmark simulation of KB contains 200,000 particles and starts up with 
222 Jeans masses, or nj^nit = 901 particles per Jeans mass. The creation of sink particles requires an increase 
in density by a factor of 40,000, corresponding to a decrease in Jeans mass by a factor 40, 000 1 / 2 = 200. 
Hence, by the time sinks arc created, the Jeans mass is down to 901/200 ~ 5 particles, which is clearly 
insufficient to resolve it. 

Scenario II also describes simulations without particle splitting, or with an insufficient number of splitting 
generations, but differs from Scenario I in the choice of the accretion radius r acc . In this scenario, r acc is set 
to a value larger than (Rj) c , by requiring that each sink particle must be made of at least nj im i n particles, 
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as described in §2.4. As a result, fewer sink particles are created, and their initial masses (that is, before 
they grow by accretion) exceed the Jeans mass by a factor Tz t / 5m in/(^j)c- 

Finally, Scenario III describes simulations with a sufficient number of splitting generations. Particle 
splitting prevents artificial fragmentation, and when sink particles are created, each Jeans-mass fragment 
contains enough particles to be replaced by a sink particle, without lumping fragments together. The number 
of sink particles formed under Scenarios I and III should be comparable, since each sink particle is created 
with a mass M ~ (Mj) c . However, in Scenario III, the Jeans mass was fully resolved throughout the entire 
calculation, thus preventing artificial fragmentation, while in Scenario I, artificial fragmentation leads to 
sub- Jeans-mass fragments, which presumably get lumped together when sinks are created. 

Notice that Scenario III is the only one we regard as satisfactory. Scenarios I and II suffer from artificial 
fragmentation. Also, in Scenario I undcrresolved fragments are turned into sinks, while in Scenario II Jeans- 
mass objects (M > Mj) are not allowed to form. 

4. THE SIMULATIONS 

4.1. Initial Conditions 

Our method for generating initial conditions is similar to the one used by KB, and is based on the 
Zel'dovich approximation commonly used for cosmological simulations. We assume that the initial density is 
described by a Gaussian random field with a density power spectrum P(k), where k is the wavenumber. The 
details are given in Appendix A. As in KB, we initially consider a power spectrum that follows a power law, 
P{k) cx k~ n , with n = 2. The rms density fluctuation S\ at scale A ~ 1/k is given by S\ ~ fc 3 / 2 P 1 / 2 (fc) ~ 
fcV 2 ~ A -1 / 2 . Hence, the density fluctuations are larger at smaller scales. 



4.2. Numerical Parameters 

For all simulations presented in this paper, we use N = 64 3 = 262, 144, n,/. m i n = nj >a .f. = 100, p c /p = 
40, 000, and / sp Ht = 8. We allow up to 2 generations of particle splitting. With A gon = 2, equation (18) 
gives Nj = 839 to be the maximum number of Jeans masses the system can contain initially. This is a 
rather large number, and it might be desirable to use a smaller one to increase the resolution per fragment, 
while retaining good statistics. We shall use instead Nj — 500 in all simulations. The initial number of 
particles per Jeans mass is then nj- m i t — N/Nj = 524. We performed three simulations, Runs A, B, and 
C, with identical initial conditions and Ag en = 0, 1, and 2, respectively, as indicated in Table 1. All other 
input parameters, including p c /p, are the same for all runs, but the accretion radius r acc is adjusted such 
that the sinks, at the time of their creation, contain at least 1 Jeans mass, and at least enough particles to 
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Fig. 2. Evolution of the number of particle per Jeans mass, nj, as the density p increases, for three 
different values of N gcn , (no splitting), 1, and 2, as indicated. The horizontal dotted line indicates nj ja .f., 
the minimum value of nj required to prevent artificial fragmentation. Particle splitting, when allowed, occurs 
when n j drops below nj >a .f., and causes nj to increase by a factor of / sp iit = 8. The vertical dotted line 
indicates the threshold density p = p c for creating sink particles. The solid dots on that line indicate the 
number of particles that end up inside each sink particle at the time of its creation, according to the various 
Scenarios (I, II, and III). 



be properly resolved (that is, n acc = max[(nj) c , nj jm j n ]). With this particular choice, Runs A and B, which 
do not have enough generations of particle splitting, will follow Scenario II; Run C, with (N gCD = 2), will 
follow Scenario III. 

Figure 2 shows the evolution of the number nj of particles per Jeans mass as the density increases. This 
number decreases as p" 1 ^ 2 , and drops below nj j& .f. = 100 when the density reaches p — (524/100) 2 = 27.5. 
Without particle splitting (N gcn = 0), nj will drop down to (nj) c ~ 3 when the density reaches p = p c . The 
accretion radius r acc is set to a value larger that (Rj) c , such that groups of nj ;ln i n = 100 particles will be 
converted into sink particles of mass M s ; n k ~ 38Mj. With iV gon = 1, particles will split once their density 
reaches p — 27.5 (first star in Fig. 2), increasing nj from 100 to 800. Then, nj will keep dropping as p 
increases, reaching (nj) c <~ 21 at p = p c . Groups of 100 particles will then be converted into sink particles 
of mass Msink ~ 4.8Mj, Finally, with TVgon = 2, particles will split a second time when the density reaches 
p = 1759 (second star in Fig. 2), increasing nj to 800 again. Eventually, nj will drop to (nj) c — 168 when the 
density reaches p = p c and groups of 168 particles will be converted into sink particles of mass M S i n k ~ Mj. 
Hence, in a simulation with a sufficient number of splitting generations, the number of particles nj in a Jeans 
mass follows a seesaw pattern, always staying above nj ja .f and thus preventing artificial fragmentation; Jeans 
masses are always properly resolved. 
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The accretion radius r acc must be determined experimentally. We ran the code with various values of 
^acc, up to the point when a few sinks (20 or so) have formed. We then check how many gas particles were 
removed when each sink was created. There is an "optimal" number of particles, which is the maximum of 
(nj) c and nj ;ln i n (100 for Runs A and B; 168 for Run C). If the number of particles exceeded significantly 
that optimal number, this indicated that r acc was too large. We would then try with a smaller value, and 
iterate until the number of particles turned into each sink was close to the optimal number. Notice that it 
could not be smaller, because the code would then delay the formation of sinks until enough particles have 
fallen inside r acc . 

Once r acc is determined, we must ensure that the resolution of the algorithm is sufficient to resolve that 
length scale. For the hydrodynamics, this is achieved by allowing the smoothing length of the SPH particles 
to shrink down to values smaller than r acc in dense regions. For the gravity, the particle- mesh part of the 
P 3 M algorithm uses a 128 3 grid to calculate the gravitational force. The corresponding length resolution 
is about 2 Ax « 0.016Lb ox , where Ax = L^ ox /128 is the cell size. The short-range correction part of the 
algorithm extends the resolution below the cell size, and the softening length r m i n can be chosen arbitrarily. 
We choose r m ; n to be slightly smaller than r acc . With these particular choices, the gravitational force will 
be properly resolved at all scales down to the scale corresponding to sink formation. The hydrodynamical 
forces will also be properly resolved, as long as the smoothing lengths shrink down to a value comparable 
to r m ; n or less. This will happen only if the mass resolution is large enough, that is, a Jeans mass contains 
at least nj ia .f.. If this is not the case, however, the hydrodynamical forces will be underresolved, and this is 
precisely what causes artificial fragmentation. 

In low density regions where the smoothing lengths h are much larger, there is clearly a mismatch 
between the resolutions of the gravity and the pressure force. This should not matter much, since particles 
are widely separated in these regions, and both gravitational and pressure forces are properly resolved at 
that scale. We have not observed any artificial clumping of SPH particles in low density regions resulting 
from the gravity being overresolved compared to the hydrodynamics (see Fig. 5 below). 

Table 1 lists the values of (nj) c , n Jymin , n acc , r acc , r min , and the number (iVj) si nk = n acc /(nj) c of Jeans 
masses inside sink particles at the time of their creation. Notice that the simulation with -/V gon = 2 is the 
only one that forms Jeans-mass objects [(iVj) s i„k = 1]. 



4.3. Computational and Physical Units 

The calculations are performed in computational units. In these units, the total mass M tot of the system, 
the box size Lbox, and the gravitational constant G are unity. Hence, the mean density p is also unity, and 
time is expressed in units of {Gp)^ 1 / 2 . The initial Jeans mass is given by 



M ^ = {2G-p) br) 



3/2 /A _ x -1/2 

(19) 



where k is the Boltzman constant, and p is the mean mass per molecule (Tohline 1982). In computational 
units, this reduces to 

M J;init = (ly'Vi-V'' = 1.0513c 3 / 2 , (20) 
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where e = 3kT/2p for a gas with 7 = 5/3. Since Mj^ n i t = M tot /Nj — 1/Nj, this reduces to 

e = 0.9672AJ 273 . (21) 

This is the prescription for choosing the initial value of e, which then remains constant throughout the 
calculation under the assumption of isothermality. With our particular choice of Nj = 500, we get e = 
0.01535. 

The Jeans criterion is given by 

> nj,a.f. • 22 

m 

We set Mj — Mj^^p/ p)^ 1 / 2 = Mj^n/p 1 / 2 , and eliminate Mj } \ n \ t using equation (20). We get 

^^125^- ™ 

Whenever the internal energy, density, and mass of a particle violate this condition, that particle is split by 
the algorithm (unless this would exceed the maximum number of splitting generations -/V gen allowed). 

The simulations are scale- free, and can be rescaled to any physical size of interest. To rescale to physical 
units, we choose particular values for the temperature T and mean density p of the cloud, and compute the 
initial Jeans mass Mj^ n - lt using equation (19). The total mass of the system is then M to t — NjMj,i n n, where 
Nj was the value used in equation (21) to set the initial conditions, the size of the box is ibox = (-Mtot/ p) 1 ^ 3 , 
and the physical time is obtained by multiplying the computational time by (Gp) -1 / 2 . This will be illustrated 
in §7 below, where we rescale the results of our simulations to several physical densities and compare to 
particular physical systems of interest. 



5. RESULTS 

Notation in this field is not standardized. For this paper, we adopt the following terminology, partly 
to facilitate comparison to the work of KB. A "dense region" is a part of a molecular cloud that is likely to 
form stars; for this paper, we focus on massive, dense regions able to form clusters. A "clump" is a region of 
enhanced density in the original dense region. A "core" is a clump that has become gravitationally unstable. 
We will later identify cores with sink particles, and we will interpret the mass function of cores in terms of 
the mass function of stars and substellar objects. Note that observers often use the work "core" to describe 
what we call here "dense regions." 



5.1. Global Evolution 

Figure 3 shows the evolution of the system for Run A (A gcn = 0). Gas particles are shown in blue, 
and sink particles in black. Following KB, we shall identify these sinks as protostellar cores, and from 
now on we shall use the term "sink" only when discussing the properties of the algorithm, rather than the 
physical interpretation of the results. Each panel is accompanied by three numbers, the time, the number 
of cores, and the percentage of mass accumulated in cores. The system evolves rapidly into a network of 
intersecting filaments. Around t = 1, cores start to form inside the dense knots located at the intersections 
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Fig. 3. — Global evolution of the system, for Run A. Each box shows a snapshot of the system, with SPH 
gas particles represented by blue dots (for clarity, only 1/8 of the gas particles are plotted), and protostellar 
cores represented by large black dots. For each snapshot, the time in units of (G/5) -1 / 2 is indicated in the 
top right corner, while the numbers in the top left corner indicate the number of cores and the fraction of 
the total gas that has been accumulated into cores, respectively. 
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Fig. 4. — Final state of the system, for Run A (left), B (middle), and C (right). The symbols and labels 
have the same meaning as in Figure 3. 

of the filaments. By t — 1.2, the densest knot already contains a large cluster of cores, and the filaments 
themselves start to fragment into cores. By t = 1.6, most of the filaments have disappeared, having turned 
into cores, and these cores are falling toward the main cluster. By t = 2, 90% of the gas has been turned into 
cores, and a dense cluster of cores has formed. During the later stages of the evolution, the remaining gas 
gets accreted onto the cores, and the cluster starts to evaporate: cores are ejected by close encounters, while 
the reminder of the cluster contracts and gets more tightly bound. All calculations terminate at t = 2.5, 
when 97%-99% of the gas has been turned into cores. 

Figure 4 shows the final state of the system, at t = 2.5, for all three runs. The final location of the cluster 
is about the same for all three runs, but the total number of cores significantly increases as we allow more 
splitting generations. The final number of cores is 213, 964, and 2876 for Runs A, B, and C, respectively. 
Also, there is less gas remaining in the system for Runs B and C, compared to Run A. 

Figure 5 shows a series of zooms at an early time slice (t = 0.9030) for Run C (iV gon = 2). The black, 
green, and blue dots represent the generation (unsplit), generation 1 (split once), and generation 2 (split 
twice) particles. The masses of these different types of particles have ratios (64:8:1). As explained in §4.2 
above, the first splitting occurs at density p = 27.5 and the second one at density p = 1759. Hence, the 
boundary between the black and green particles is an isosurface of density p = 27.5, and the boundary 
between green and blue particles is an isosurface of density p — 1759. Notice that these boundaries are 
quite sharply delimited, indicating that nothing peculiar is happening there. Indeed, the density varies 
smoothly across this entire region, and the transition between particles of widely different masses produces 
no observable feature in the density profile. 

The bottom right panel of Figure 5 zooms in on the intersection of filaments to show a concentration 
of sinks (red circles), which have formed earlier and have moved toward each other as the density increases. 
In this region, a gas particle might be located inside the accretion radius of several sinks, in which case it 
will accrete onto the sink for which the total energy is the lowest (that is, the binding energy is highest), 
as explained in §2.4. The bottom left panel shows a further zoom-in on the filament to the lower left. In 
this particular case, the cylindrical collapse of a filament leads to the formation of a chain of sink particles 
separated by a distance equal to twice the accretion radius r acc , reminiscent of observations of dense regions 
within filaments (Schneider & Elmegreen 1979). The circle shows the accretion radius for a particular sink. 
Gas particles that enter that accretion radius will be accreted by the sink (most particles appearing inside 
the circle in Fig. 5 are seen in projection). 
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Fig. 5. — Early time slice [t = 0.903) of the evolution of the system, for Run C. Top left panel: the entire 
computational volume. For clarity, only 1/8 of the particles are shown. Top right panel: zoom- in on a dense 
region at the intersection of two emerging filaments (with all particles shown). The colors black, green, 
and blue corresponds to particles of generation 0, 1, and 2, respectively. Sink particles are shown in red. 
Bottom right panel: zoom-in on the central, high-density region. Cavities around sink particles are visible. 
Bottom left panel: Zoom-in on a filament. The cavity around one sink particle is shown by a circle of radius 
^acc = 0.0005. Particles appearing inside the cavity are either background or foreground particles. 

5.2. Particle Splitting and Sink Particles 

Figure 6 shows stacked histograms of the number of particles within each generation (left panels) and 
the total mass contained within each generation, and the division of mass between particles and sink particles 
(right panels). In the absence of particle splitting (N gcn = 0, Run A, top panels), the number of particles 
remains constant until sink particles are created, and then steadily decreases as sink formation and accretion 
onto sinks removes particles from the simulation. When particle splitting is included (Runs B and C, middle 
and bottom panels), the total number of particles AT to t initially increases as particles split, but this process 
competes with accretion onto sinks, and eventually N to t decreases. In Table 2, we list, for each run, the initial 
number of particles N- U1 i t , the maximum number of particles iV pca k reached during the simulation, and the 
effective number of particles N e g, defined as the number of particles a simulation without particle splitting 
would need to have the same resolution (N e {[ = Mnit/j^ut")- Runs A, B, and C have effective resolutions of 
64 3 , 128 3 , and 256 3 particles, respectively. Of particular interest is the ratio N pea k/N s, given in the last 
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Fig. 6. — Left panels: stacked histograms of the number of particles of various generations, versus time. 
The labels 0, 1, and 2 indicate the various generations of particles. The top curve on each panel shows the 
total number of particles, N to t- Right panels: stacked histograms of the mass contained in each generation 
of particles, and in sinks, versus time. Top panels: Run A; middle panels: Run B; bottom panels: Run C. 



column of Table 2. This ratio measures the "efficiency" of the particle splitting algorithm. As the number 
of splitting generations increases, the peak number of particles becomes significantly lower than the effective 
number of particles, resulting in a substantial saving of computational effort relative to a simulation without 
particle splitting. From Table 2, we can infer that every additional splitting generation makes the algorithm 
more "efficient" by a factor of order 4 — 5. 

At this point, we need to discuss a possible alternative to particle splitting, called the "zoom-in" ap- 
proach, which is often used in numerical simulations. This approach would consist of first running a low- 
resolution, N = 64 3 simulation, identifying in the final state of the simulation the particles located in "regions 
of interest," where higher resolution would be desirable, going back to the initial conditions, replacing these 
particles only by a set of smaller particles, and then redoing the simulation. This approach would fail in 
the present case for the following reason: As we shall see in the next section, about half of the gas par- 
ticles in the simulation are eventually converted into sinks, while the other half get accreted onto existing 
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sinks. While high resolution might not be necessary for the particles that are accreted onto sinks, it is cer- 
tainly necessary for the particles that are converted into sinks, in order to prevent artificial fragmentation. 
This means that, in the initial conditions, one half of the particles would have to be replaced by a cube of 
4 x 4 x 4 = 64 particles in order to provide sufficient resolution. The total number of particles would then 
be N = 64 3 [l/2 + (4 3 — l)/2] = 8, 388, 608, or 256 3 /2. Hence, using a zoom-in approach would only reduce 
the number of particles by a factor of 2 compared with the brute force approach. The zoom-in approach is 
useful in situations where the "regions of interest" contain a small fraction of the total mass of the system. 
Here, these regions contain 50% of the mass, making this approach inefficient. 

5.3. Regimes 

Figure 7 shows the time-evolution of the number of cores and the total mass accumulated in cores. 
By comparing the time-evolution of these two quantities, we can identify four distinct phases during the 
evolution of the cloud, with each phase corresponding to a different dynamical regime. 

The initial phase corresponds to the growth regime. The initial density fluctuations grow by gravitational 
instability, to form a network of dense, intersecting filaments. This phase of the evolution terminates when 
the first cores form. Since the density threshold p c for sink formation is chosen arbitrarily, the end of this 
initial phase is also arbitrary. However, in the absence of sinks (and with infinite resolution), fragments 
would collapse and reach infinite densities at a finite time that would exceed the time of sink formation only 
slightly. Hence, choosing the time of formation of the first cores as corresponding to the end of the growth 
phase is not an unreasonable choice. 

After the first cores form, the evolution of the cloud enters the second phase, which corresponds to the 
collapse regime. To understand this regime, and the following one, we must consider a fundamental property 
of Gaussian initial conditions: the filling factors of underdense (p < p) and overdense (p > p) regions are 
initially equal, both being 1/2 of the total computational volume. Since the initial density is nearly uniform 
(p w p), about 1/2 of the gas starts up in overdense regions, while the other 1/2 starts up in underdense 
regions. In the collapse regime, the gas that started up in overdense regions collapses and is converted into 
cores. This regime is characterized by an increase in both the number of cores and the mass inside cores, 
at rates that are roughly proportional. This phase terminates when 50% of the gas, essentially the gas that 
started up in overdense regions, has been turned into cores. 

The evolution of the cloud then enters the next phase, which corresponds to the accretion regime. The 
uncollapsed gas remaining in the cloud is the gas that started up in underdense regions. The main tendency 
of this gas is not to collapse onto itself and form new cores, but rather to accrete onto the existing cores. 
As a result, the mass in cores keeps increasing, while the number of cores levels off. During this phase, the 
mass in cores nearly doubles, while the number of sinks increases by less than 30%. 



Table 2. Number of Particles. 



Run 


iVinit 


N cB 


Npeak 


Afpoak/iVoff 


A 


262,144 


262,144 


262,144 


1.000 


B 


262,144 


2,097,152 


467,636 


0.223 


C 


262,144 


16,777,216 


1,067,610 


0.064 
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Fig. 7. — Time-evolution of the mass accumulated in protostellar cores (dashed curve, left axis) and the 
number of cores (solid curve, right axis) versus time. The dotted lines separate the various phases of the 
evolution, with the corresponding regimes labelled on top. 

Finally, once most of the gas (90% or so) has been accumulated in cores, the evolution of the cloud 
enters the fourth and last phase, which corresponds to the N-body regime. The hydrodynamics becomes 
irrelevant, and the evolution of the system is governed by gravitational many-body dynamics. 

5.4. Formation and Growth of Protostellar Cores; Low-Resolution Simulation 

Figure 8 shows the final masses of the protostellar cores as a function of their birth rank. We define the 
birth rank such that the n th core formed in the simulation has a birth rank of n (KB use the term "index"). 
The top panel shows the results of Run A, without particle splitting, which can be directly compared to the 
results of KB. As KB noted, the mass tends to decrease with increasing birth rank, simply because cores 
formed later have less time to grow by accretion. However, the trend we observe is much less pronounced 
than the one found by KB. The mass range is comparable, and the most massive cores all formed early 
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Fig. 8. — Mass of protostellar cores versus birthrank. (a, top panel) Run A (N gcn — 0); (b, middle panel) 
Run B (-/V gon = 1); (c, bottom panel) Run C (iVg Cn = 2). The various symbols (crosses, asterisks, and solid 
circles) identify particular cores that are discussed in the text. 



(with the most massive one being the very first one that formed), but the distribution shown in Figure 8a 
reveals several low- mass cores that formed early, as well as several high- mass "peaks" with high birth ranks. 
Another surprising result, for Run A, is the final mass of the most massive core, which is 35% of the total 
mass of the system. This core, the first one to form, clearly experienced runaway accretion. By contrast, 
the most massive core formed in any simulation of KB has a mass of order 6 — 8% of the total mass of the 
system. It is important to understand these various features. For the rest of §5.4, we will discuss the results 
from Run A; Runs B and C will be discussed in §5.5. 
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5.4-1- Early Ejections 

The presence of cores in the bottom left corner of Figure 8a is easily explained. Some cores that form 
early can be ejected after experiencing close encounters with other cores. Once a core is ejected from the 
dense structure where it originally formed, it finds itself in a low-density region where there is little gas to 
accrete. The mass growth of that core then levels off at a constant value, a process discussed by KB. To 
illustrate this, we focus on 3 early-forming cores that end up with very low masses, cores #8, #10, and #14 
(crosses in Fig. 8a). Figure 9a shows the distribution of cores in the system at t = 1.084, just before core #10 
is ejected. The system at that time contains 19 cores, 10 of them forming a dense cluster embedded in a 
common gas filament (not shown). This cluster, indicated by the small square, is enlarged and displayed in 
Figure 9b. Figure 9c shows the trajectories of these 10 cores, from the locations where they formed (open 
circles) to their locations at t = 1.154 (solid circles). Several cores are ejected, including the cores #8, #10, 
and #14. 

Figure 9d shows the mass evolution of these 3 cores. After formation, their masses grow rapidly by 
accretion, until they are ejected. After being ejected, cores #10 and #14 no longer accrete gas, and remain 
at constant mass throughout the reminder of the simulation. Core #10 is ejected at very large velocity (see 
the quasi-straight trajectory in Fig. 9c), and moves only through low-density regions thereafter. Core #14 
remains bound to the main cluster, but ends up orbiting the cluster at a "safe" distance, never coming close 
enough to accrete gas from the cluster's envelope. Core #8 is ejected at low velocity, turns back, and falls 
into a dense clump at t = 1.541. Its mass then increases slightly by accretion. At t = 1.904, that clump, 
which still contains core #8, merges with the main cluster. The mass of core #8 increases again by accreting 
gas from the main cluster, before being ejected a second time, after which its mass remains constant. 

5-4-2. Local Competitive Accretion 

To understand the origin of the massive cores at high birth rank, we need to consider the nature of 
competitive accretion between cores. This process was described in detail by Bonnell et al. (1997, 2001). 
There are 4 basic arguments for why cores that form earlier should reach higher masses: (1) since they form 
earlier, they have more time to grow by accretion, (2) the cores that form early will deplete their surroundings 
by accretion, reducing the amount of gas available to cores that form later, (3) by being more massive, the 
cores that form early have a stronger gravitational potential than the ones forming late, making them more 
efficient in accreting the reminder of the surrounding gas, and (4) if several cores of different masses are 
present inside a gaseous clump, the most massive cores will tend to reside in the center of the clump where 
there is presumably more gas to accrete. While argument (1) is general, arguments (2), (3), and (4) are 
valid only if cores form out of the same clump, and are therefore competing for the same surrounding gas. 
If the final cluster forms by the merging of dense clumps, and if cores form in these clumps prior to the final 
merging, then there will be competitive accretion within each clump, but not across clumps. The very first 
core that forms in a particular clump can then grow by accretion and reach a high mass, no matter how late 
core formation in that clump started. 

To illustrate that, we focus on 3 cores, cores #141, #151, and #156, which are located in two late 
"peaks" (asterisks in Fig. 8a). The top panels in Figure 10 shows the system at t = 1.235, immediately 
after the formation of core #141. As we see, that core did not form inside the main cluster, but inside an 
emerging filament located away from the main cluster, and was the very first core to form there (core #5 is 
a fast-moving core that was ejected early and is seen in projection). The bottom panels in Figure 10 show 



- 26 - 



t = 1.084 



t = 1.084 




t = 1.154 




time 



Fig. 9. — (a, top left panel) snapshot of the distribution of cores in the computational volume at t = 1.084. 
The system contains a total of 19 cores; (b, middle left panel) enlargement of the region indicated by a 
square in the previous panels, showing a dense clusters composed of 10 cores; (c, top right panel) the same 
cluster is shown at t = 1.154 (solid dots), along with the trajectories of the cores between t = and t = 1.154 
(curves). The cores #8, #10, and #14, which are ejected, are indicated; (d, bottom panel) time evolution 
of the mass of cores #8, #10, and #14. 



the system at t = 1.265, immediately after the formation of core #156. Again, cores #151 and #156 formed 
in a filament away from the main cluster, and were the first two cores to form there. Local competitive 
accretion enables these cores to grow and reach a high mass before they fall into the main cluster. 



5.4-3. Violent Infall and Late Starburst 

There is another peak in Figure 8a that we wish to consider. This peak contains cores with birth ranks 
199-203 (solid circles in Fig. 8a), that formed very late. These 5 cores formed almost simultaneously, between 
t = 1.595 and t = 1.609, all inside a small region of diameter 0.006 located right in the middle of the cluster. 
Figure 11 shows a histogram of the formation time of cores. The formation of cores #199-203 in Run A 
appears as a burst, which is indicated by an arrow in the top panel. We need to understand how this burst 
occurred and how these cores managed to reach a high mass in a such a crowded environment. Clearly, 
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Fig. 10. — (a, top left panel) distribution of cores at t = 1.235, immediately after core #141 formed; (b, 
top right panel) enlargement of the region indicated by a square in the previous panel, showing both gas 
particles (small dots) and cores (large dots). Cores #5 and 141 are indicated by numbers; (c, bottom left 
panel) and (d, bottom right panel) similar to (a) and (b), but at time t — 1.265, immediately after core 
#156 formed. 

we cannot invoke local competitive accretion here, since these cores, unlike cores #141, #151, and #156, 
formed in a region that already contained many other cores that were already significantly more massive. 

Consider the state of the system at time t ~ 1.5 — 1.6 (see 8 th and 9 th panels in Fig. 3). At that stage, 
there is still a substantial amount of gas located far from the main cluster. Some of that gas, which is located 
in dense filaments, will form new cores, that will later fall into the main cluster. The remainder of the gas, 
which is located between the filaments, does not reach high enough densities to trigger the formation of 
cores. It will therefore remain in the form of gas until it falls into the main cluster. Because that gas comes 
from large distances, it gets accelerated over a long period of time, and therefore falls inside the main cluster 
at large velocities. This large velocity reduces the effectiveness of gravitational focusing, making it more 
difficult for that gas to accrete onto the cores already present in the cluster. Instead, that gas is rapidly 
decelerated by a strong shock and gets compressed to a very high density, which triggers the formation of 
several cores in a burst. 
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Fig. 11. — Histogram of the number of protostellar cores versus their epoch of formation, (a, top panel) Run 
A (iV gon = 0); (b, middle panel) Run B (N gcn = 1); (c, bottom panel) Run C (iVg 0n = 2). The arrow in the 
top panel identifies the late starburst. The first three cores formed in Run A are also identified. 



5.4-4- Runaway Accretion 

To understand the very large mass of core #1, we need to consider the formation history of the first 
few cores. Looking again at Figure 11a, we immediately notice something special about the first three cores. 
These cores formed at well-separated times, after which core formation proceeded very rapidly. Table 3 
shows, for the first 5 cores, the quantity irii(tj), defined as the mass of core i at the time of formation of 
core j (with j > i). All cores are created roughly with the same mass (diagonal in Table 3). However, 
core #2 formed significantly later than core #1, and during the time interval ti—t\, core #1 accreted more 
than 20 times its original mass (growing from 0.000435 to 0.009205), so that when core #2 formed, core #1 
was already more massive by a factor of 23. The process then repeats itself: core #3 formed significantly 
later than core #2, and during the time interval <3 — £2, core #1 and #2 experienced significant growth. By 
t = t 3 , these 3 cores have mass ratios 55:10:1. This process goes on, and by t = i 4 , the first 4 cores have mass 
ratios 77:11:10:1 (core #2 is ejected between t 3 and £4, and stops growing afterward). Finally, this process 
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terminates: core #5 is created almost immediately after core #4, providing little time for the latter to grow, 
so that at t = t 5 , cores #4 and #5 have comparable masses. We are then in the "big peak" in Figure 11a, 
with many cores forming at comparable times. 



5.5. Formation and Growth of Protostellar Cores; Higher-Resolution Simulations 

We have focused so far on Run A, and investigated the origin of the various features found in Figure 8a. 
Figures 8b and 8c show the results for Runs B and C, respectively. The results are qualitatively similar to 
Run A. There is a definite tendency of the mass to decrease with increasing birth rank. The distributions 
are very wide, however, and show several peaks at high birth rank, which we again attribute to competitive 
accretion. 



5.5.1. Early Ejections 

There is a noticeable difference in the higher-resolution runs: the absence of very-low mass cores at small 
birth ranks. Early ejections seem far less common in these simulations. This can be explained using the 
following argument. Consider an envelope of gas of mass M and radius R, containing several cores of mass 
to, forming a small, bound cluster of radius r, inside which the mean core spacing is A. The gravitational 
force / between cores is proportional to m 2 /A 2 , and therefore the acceleration a of the cores varies with to 
as 

f m 2 /A 2 to , . 

a = J- cx —L cx — . 24 

TO TO A^ 

For a fixed cluster radius r, the core spacing A depends on the number of cores, and for a fixed cluster mass, 
that number of cores depends on the resolution, such that A cx to 1 / 3 (if the mass per core is smaller, the 
cluster contains more cores, thus those cores are closer to one another). Hence a cx to 1 / 3 . If a core is ejected 
from the cluster, its terminal velocity should be of order v ~ (2ar) 1 / 2 cx to 1 / 6 . But the escape velocity from 
the whole system is of order v csc « (GM/i?) 1 / 2 , independent of to. Therefore, lowering to makes it more 
difficult for cores ejected from the cluster to escape the common envelope of gas. 



5.5.2. Local Competitive Accretion and Late Starburst 

Figures 8b and 8c show several peaks, corresponding to high- mass cores with high birth ranks. We 
examined several cases, indicated by asterisks. In all cases, we found that the high masses resulted from 

Table 3. Mass of the First Cores (xlO 3 ). 



Mass 


ti = 0.9010 


*2 = 0.9522 


t 3 = 1.0055 


U = 1.0356 


t 8 = 1.0361 


mi 


0.435 


9.205 


21.942 


31.693 


31.925 


m.2 




0.405 


3.922 


4.494 


4.494 


TO3 






0.401 


4.166 


4.242 


m.4 








0.412 


0.435 


TO 5 










0.381 
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local competitive accretion. All these cores formed in regions located away from the main cluster, often in 
emerging filaments. In several cases, these filaments contained a string of cores (as seen, for instance, in 
Fig. 5), and the core that reached a high mass was the one located at the very end of the string. This is an 
example of local competitive accretion, and perhaps the core at the end is able to accrete from a larger solid 
angle. 

We did not find a late starburst (or coreburst) in the high-resolution simulations. The gas located in 
very low regions does fall into the cluster in the form of gas, but tends to be accreted onto the cores rather 
than forming new cores. Accretion is much more efficient in the higher-resolution simulations simply because 
of the sheer number of cores in the cluster at late time. That number is very small for Run A because (1) 
core #1 contains more than half the final mass of the cluster, and (2) as we shall see in §5.6, the cluster suffers 
a great deal of evaporation at late time in Run A. Hence, the cross section for accretion onto cluster cores 
is greatly reduced in the low-resolution simulation compared with the high-resolution ones, which explains 
why the late starbust phenomenon is seen only at low resolution. 

5.5.3. Runaway Accretion 

We did not observe the kind of runaway accretion experienced by the first core in Run A. While that 
core reaches a final mass of 0.358 (35.8% of the total mass in the computational volume), the final masses of 
the first cores formed in Runs B and C are 0.00295 and 0.000803, respectively. Actually, the core which ends 
up with the largest mass in Run C is not core #1, but rather core #707, with a final mass of 0.00263. In 
Run A, runaway accretion occurred because the first core formed significantly earlier than the others. With 
higher resolution, several lower-mass cores form almost simultaneously inside the first region that reaches 
the threshold density, and competitive accretion among these cores prevents runaway accretion. 

Comparing all these results with the ones presented in §5.4, we conclude that local competitive accretion 
is a fundamental process that occurs at all resolutions, except in situations where a small number of cores 
results in the formation of a single cluster (as in the simulations of KB). In such a case, competitive accretion 
does occur, but it is not local. Other phenomena, like early ejection, late starburst, and runaway accretion, 
appear to be peculiarities of low-resolution runs, and become less common as the resolution increases. 

5.6. The Protostellar Cluster 

All three simulations end up forming a dense cluster of protostellar cores. In this section, we study the 
assembly history and structure of that cluster. 

5.6.1. Mass and Cluster Members 

To identify the cores that belong to the cluster, we first use visual inspection to make an initial estimate 
(r c i)init of the center of the cluster. Then, using an iterative method, we find a self-consistent solution for 
the center of mass r c i and radius r 2 oo of the cluster, such that (I) the center of mass of the cores located 
inside a sphere of radius r 2 oo centered on r c i is indeed r c i, and (2) the mean density P200 inside that sphere 
is 200 times the mean density in the system. Notice that these densities are computed using the cores only, 
without the gas, as in KB. 
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Fig. 12. — (a, top panel) evolution of the number of cores inside the main cluster; (b, bottom panel) evolution 
of the mass in cores inside the main cluster. On both panels, the three curves correspond to Run A (N gen = 0, 
solid curve), Run B (N gCD = 1, dotted curve), and Run C {N gCD = 2, dashed curve). On the top panel, the 
number of cores for Runs A and B have been multiplied by factors of 25 and 2.75, respectively, to allow 
a better comparison. The thick vertical lines indicate the transitions between the collapse, accretion, and 
N-body regimes, as indicated. 



Figure 12a shows the time-evolution of the number of cores in the cluster, for all three runs. We have 
rescaled the number of cores for Runs A and B by factors of 25 and 2.75, respectively, for comparison 
with Run C. The number of cores for Run A (solid curve) varies tremendously, in a somewhat chaotic way. 
These large fluctuations can be attributed to small statistics, Run A having the smallest total number of 
cores. However, this is not the only explanation. Our method for determining cluster membership assumes 
spherical symmetry, but in Run A, the cluster tends to be triaxial during most of the calculation, and also 
the sphere of radius r2oo which contains the cluster members tends to follow the motion of core #1, since 
that core contains most of the mass of the cluster. So as core #1 experiences brownian motion around the 
center of the cluster, cores located near the surface of the sphere keep "falling" in and out of the cluster. 
This does not occur with the higher-resolution runs, because the cluster tends to be quite spherical, and no 
single core dominates its mass. 

During the collapse regime (t < 1.34), cores are forming inside the cluster, increasing the number of 
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cluster members. There are few ejections during this phase, and the net effect is a near-monotonic increase 
in the number of cores. Once the system enters the accretion regime, few new cores are forming, and the 
competing processes are ejections from the cluster and accretion of cores that formed in secondary clumps 
outside the main cluster. Indeed, the main cluster experiences a major merger with another, comparable 
cluster, at t ~ 1.29, near the end of the collapse phase, resulting in a sudden increase in the number of 
cores. For Runs B and C, the number of cores tends to increase during the accretion regime, indicating that 
accretion of cores formed in subclumps dominates over ejections. For Run A, however, the ejections tend to 
dominate. Once the system enters the N-body regime (t > 2.02), the number of cores steadily decreases for 
Runs A and B. All the cores formed outside the main cluster have been accreted along with the remaining 
gas, and ejections become the dominant process. However, the number of cores in Run C remains nearly 
constant, indicating that very few ejections are taking place. 

Figure 12b shows the evolution of the mass of the cluster. The results are nearly identical for all three 
runs up to t ~ 1.29, when the major merger takes place. From that point, Run A differs significantly from the 
other two runs, as the mass of the cluster grows much more slowly. Interestingly, in the N-body regime, the 
mass of the cluster, for Run A, remains nearly constant even though the number of cores drops significantly. 
This indicates that only low-mass cores are ejected at late time. 

The most striking result shown in Figure 12b is the similarity between the results of Runs B and C. 
The curves are essentially indistinguishable up to the beginning of the accretion phase, and then remain 
very similar up to the end of the simulation, having in particular the same local maxima at t — 1.33, 1.58, 
and 2.05, and the same local minima at t — 1.37 and 1.61. Notice also the sudden increase in mass at 
t = 1.54, corresponding to a merger with a smaller subcluster of cores. This merger does not occur in the 
low-resolution simulation. 



5.6.2. The Density Profile 

We computed the density profile of the cluster by adding up the mass of the cores inside radial bins. 
Figure 13 shows the evolution of the profile, for all three runs. For Run A, the profile is very steep at 
the center, simply because core #1, which undergoes runaway accretion, totally dominates the mass of the 
cluster. For Runs B and C, the cluster profile starts up roughly as a power law, but soon acquires a core/halo 
structure, as in the simulations of KB. In these two runs, the central density drops and the profile flattens at 
late time (t > 2), corresponding to the epoch where the total mass of the cluster stops growing (see Fig. 12). 
At these late times, the central part of the cluster inflates, probably as a result of gravitational heating by 
tight binary cores. 

Figure 14 shows a comparison among the three runs of the cluster profile at some particular times 
(t = 1.13, 1.5, and 2). The density profiles found in Runs B and C are very similar, at all epochs. This is 
quite striking, considering that the number of cores in the cluster roughly quadruples between t = 1.13 and 
t = 2, while the mass in the cluster increases by a factor of 9. The profile for Run A is totally different. 
Most of the mass of the cluster in contained in core #1, located in the center of the cluster. 
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Fig. 13. — Density profile of the main cluster, for Runs A (top panel), B (middle panel), and C (bottom 
panel). The various curves correspond to different times, from bottom to top: t = 1.00 (bottom curves), 1.13, 
1.50, 2.00, and 2.50. There is no curve plotted for t = 1.00 in the top panel because the cluster contained 
only two cores at that time in Run A. 



5.7. Initial Mass Function of Protostellar Cores 



Figure 15 shows the final mass distribution of protostellar cores at t = 2.5, when nearly all the gas has 
been accreted. For Run A, the existence of a core containing 35% of the total mass of the system results in a 
skewed distribution. Interestingly, the simulations of Tillcy & Pudritz (2004) show similar features (Fig. 17 
in their paper) , suggesting that runaway accretion is happening in their simulations as well. For Runs B and 
C, the distributions are roughly log-normal, with some skewness toward high masses, as the dashed curves 
show. KB reported that the mass distribution peaks at the initial Jeans mass of the system. Our results are 
consistent with that claim, but only for Run A. As the number of splitting generations, or equivalently the 
mass resolution of the algorithm, increases, the distribution shifts to lower masses while keeping the same 
width in log M. Clearly, it is the resolution, and not the initial Jeans mass, that determines the location of 
the peak. As the resolution increases, the algorithm can form cores with smaller masses, until the resolution 
is so high that the Jeans mass (Mj) c at the threshold density p — p c is resolved. Hence, increasing the 
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Fig. 14. — Density profile of the main cluster at times t = 1.13, 1.50, and 2.00, as indicated. The solid, 
dotted, and dashed curves correspond to Runs A, B, and C, respectively. The numbers in the bottom left 
corner of each panel indicate the number of cores in the cluster at these times (top number: Run C; middle 
number: Run B; bottom number: Run A). 

resolution moves the low-mass end of the distribution to lower masses. This leads to the formation of a 
larger number of cores, which then have to compete for accretion. As a result, the high-mass end of the 
distribution also moves to lower masses. Note that the shift to lower masses results only from increasing 
iVggn because p c is the same in all three runs. 

These log-normal distributions (for Runs B and C) are consistent with the numerical results of KB. 
However, they are inconsistent with observations that show a peak at small masses and a power-law like 
behavior at the high-mass end. There are several possible explanations for this. Once sinks form, they 
can no longer fragment. In the real world, these objects might fragment, increasing the number of low- 
mass objects relative to high-mass ones (notice that sink merging, which we also ignore, would have the 
opposite effect). We did not take turbulence into account in these simulations. Turbulence could affect 
the dynamical evolution of the system differently at different scales, affecting the final shape of the mass 
distribution. Probably even more importantly, feedback effects, if included, could slow down the growth of 
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Fig. 15. — Mass distribution of protostellar cores, (a, top panel) Run A (N gcn — 0); (b, middle panel) Run B 
{N gen — 1); (c, bottom panel) Run C {N gcn = 2). On each panel, the left and right dotted lines indicate the 
Jeans mass at densities p — p c and p — p, respectively. The dashed curves show the results of a least-square 
fit to a log-normal distribution. 



massive cores by accretion. Finally, we could consider changing the slope of the power spectrum of initial 
density fluctuation, though KB did this (Klessen & Burkert 2001), and found that the distributions of core 
masses remained log-normal. 



6. CONVERGENCE 

The three simulations we have performed all start with 64 3 particles, and use identical initial conditions. 
However, particle splitting increases the effective resolution of the simulations, which is 128 3 for Run B 
and 256 3 for Run C. Since the initial conditions are identical but the effective resolution varies, this set of 
simulations constitutes a convergence study, which can be used to estimate the minimum resolution necessary 
to obtain reliable results. 
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It is clear, however, that some results simply will not converge. By requiring that sink particles are 
created with enough gas particles to insure that they are resolved (that is, Scenario II, by opposition to 
Scenario I), the initial mass of sinks depends on resolution, and as a result the IMF of cores shown in 
Figure 15 shifts to lower masses as the resolution increases. However, the results of Run C have converged, 
since that simulation does resolve the Jeans mass. If we added a Run D with i\r gcn = 3 to our set of 
simulations, the results of Runs C and D would be identical, because a third level of splitting would never 
occur: before the density gets large enough to make particles split a third time, these particles would turn 
into sinks (see Fig. 2). This convergence is numerical in the sense that numerical parameters like p c and 
iVgcn determine the solution that the simulation converges to. A solution that would converge physically 
does not exist, and therefore cannot be achieved, in a system with an isothermal equation of state, because 
no physical process limits the minimum mass of cores. In the real universe, the assumption of isothermality 
breaks down at high densities when the gas becomes optically thick, and that in turns leads to a physical 
minimum mass for cores. 

Looking at the macroscopic properties of the final cluster of cores, it is clear that the results of Runs 
A and B are significantly different, while the results of Runs B and C are very similar, indicating that 
convergence has been achieved. In particular, the mass history of the cluster (bottom panel of Fig. 12) and 
the density profile at various times (Fig. 14) are strikingly similar for Runs B and C. In Run A, the first core 
formed underwent runaway accretion, which affected the further evolution of the cluster. No such runaway 
accretion occurred in Runs B and C. We believe that the likelihood of such occurrence is reduced as the 
resolution increases, because as more cores are formed, the time interval between the formation of a core 
and the next one is reduced, thus increasing the competition for accretion. 

Finally, there are other results, arguably less interesting, that show convergence. In particular, the 
formation time histograms shown in Figure 11 are quite similar for Runs B and C, and different for Run A, 
with the bulk of the cores forming at later time. 

7. ILLUSTRATIVE EXAMPLES 

To compare the results to observations, it is convenient to use the relations in §4.3 to convert the 
densities, masses, etc. to physical units. Because the simulations contain 500 Jeans masses initially, the 
total mass and other properties arc uniquely fixed by a choice of temperature and density. We set the 
temperature at 10 K for all these examples because this temperature is characteristic of both dust and gas 
temperatures in well shielded regions before stars form (Leung 1975; Evans et al. 2001; Young et al. 2004). 
We will consider below constraints imposed by the assumption of isothermality. Observers commonly use 
total particle density [n = n(H2) + n(Hc)], where p = p, n mHn, with itih the atomic mass unit. For a fully 
molecular cloud with 25% helium by mass, p n = 2.29. We will use n for the initial density. 

With equation (19) for the Jeans mass, we have in physical units: Mj^n = 6.33T 1,5 n _o ' 5 M0 = 
200n-°- 5 M Q for T = 10 K. It follows that M tot = 1.0 x lO 5 n"°- 5 M , and the size of the region, L box = 
(M t ot//5) 1/3 = 121n-°- 5 pc The dynamical time, t dyn = (Gp)" 1/2 = 6.3 x 10 7 n"°- 5 yr. Values for these 
quantities are given in Table 4 for different values of n. 

The assumption that the gas remains isothermal depends on its ability to cool. In dense regions, the gas 
cools by collisions with dust grains, which radiate in a continuum (Goldreich & Kwan 1974; Doty & Neufeld 
1997; Young et al. 2004). To remain isothermal, the optical depth in the continuum near the peak emission 
wavelength should be less than unity, measured from the center of the region to the edge. The initial optical 
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depth is computed from r = Kpibox/2, where k is the opacity of dust per gram of gas. Emission from dust 
at 10 K peaks at a wavelength around 350 ^im. Calculations of dust opacities for dust that has coagulated 
and acquired ice mantles, as may be expected in dense regions, have been done by Ossenkopf & Henning 
(1994). Observations are generally well matched by the opacities from column 5 in their table, known as 
OH5 opacities. The value for 350 /im, assuming a gas to dust ratio of 100, is 0.1 cm 2 g _1 of gas. The values 
of r in Table 4 are computed from these assumptions. 

Because we form sinks at a density / s i nk = 4x 10 times the initial density (using the notation of §3.1), 
we must check that the region around the sink is optically thin just before sink formation. A convenient 
measure for this is the optical depth calculated for the radius of the SPH kernel at p c = f sm kP- The 
radius of the kernel is about 3 times the local particle spacing Ar. The particle spacing is constrained 
by Ar < 0.5^/7^ n part2 Nson , where n par t is the initial number of particles along the edge of the 
computational volume (64 for these simulations). Thus, we have for the center to edge optical depth of a 
kernel, r kern = tF, where F = &fHX n vlvt 2 ~ Nsm ■ For ^gon = 0, 1, 2, F = 109.6, 54.8, 27.4. Particle splitting 
helps to keep the optical depth in a kernel from rising too far above the initial value, since increasing N gen 
lowers F. Note that Tk or n oc r oc n 5 because of the scalings compelled by assuming that we have 500 Jeans 
masses. With the criterion that Tk orn < 1, these calculations would be valid up to n ~ 10 5 cm~ 3 , so we do 
not show any entries in Table 4 with n above that value. With / s i n k = 4 x 10 4 , the maximum density would 
be about 4 x 10 9 cm~ 3 . Alternatively, we could use the accretion radius (r acc , see §2.4) instead of the kernel 
radius. These two radii turn out to be very similar, so the results are not significantly affected. 

The 4 cases listed in Table 4 are simply scaled to different values of the initial density. The first is 
roughly similar to the example given by Klessen & Burkert (2000), which they compare to the Taurus cloud. 
Our case 1 has about twice the mass of that of KB because we have roughly twice as many Jeans masses. 
Our definition of Jeans mass is slightly different, and other differences combine to make our value for Ly, ox 
about twice their value. The cases with higher initial densities approximate the conditions in cluster forming 
regions. Case 3 is similar to the conditions in the L1688 (Ophiuchus) cluster, which has 500 — 1000M Q in a 
region extending about 1 pc (Johnstone et al. 2000). Case 4 approximates conditions in the massive, dense 
regions studied by Shirley et al. (2003a); the median radius of gas traced by CS emission in their study is 
0.37 pc, and the mean virial mass is about 1200 M Q . Mass estimates based on the 350 /zm emission (Mueller 
et al. 2002) for a subsamplc of those sources yield mean values around 300 M Q . 

Putting the mass functions shown in Figure 15 into physical units is probably premature because we 
have not yet included many effects, such as feedback from early star formation. If we proceed with these 
caveats in mind, we find that the mass function, in the simulation with A gcn = 2, and for Case 3, peaks 
around 0.3 M , with a maximum mass around 2.0 M Q and a minimum mass of 0.03 M Q , well into the 
region of brown dwarfs. Case 4 produces even lower mass objects, with a peak in the distribution around 0.1 
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Mq. This is not surprising since there is no limit to how small the Jeans mass may get with an isothermal 
equation of state. Clearly the other effects, which are the subject of future work, will need to be included 
before these distributions can be compared to the observed mass functions. 

8. DISCUSSION 

We have performed three simulations of fragmenting molecular clouds, using identical initial conditions, 
but different levels of particle splitting. Each additional level corresponds to an increase by a factor of 2 
in length resolution and 8 in mass resolution. We have discovered several interesting phenomena, which we 
described in §5 above. Some of the results, such as the runaway accretion and the late starburst, were found 
only in the lowest-resolution simulation, and thus are clearly resolution-dependent. Such phenomena are 
peculiarities that happen sometimes when the number of cores is relatively small. We found three results 
that are not resolution-dependent: (1) the existence of four distinct regimes in the evolution of the cloud, 
(2) the phenomenon of local competitive accretion, and (3) the tight relationship between the mass range of 
the IMF and the resolution limit of the algorithm, until the calculation is fully resolved. 

In these initial simulations, we have neglected several physical processes, such as turbulence, non- 
isothermality, and feedback, which will be addressed in future work. We believe that the three main results 
stated above are robust, and will remain valid once we include additional physical processes (though this 
remains to be proven). First, the existence of the four distinct regimes should not be affected by additional 
physics. The gas will always end up into cores given enough time, so there will always be a N-body regime 
at the late stages. The fact that half the gas turns into cores during the growth regime, while the other 
half is accreted during the accretion regime, is most likely a consequence of the Gaussianity of the initial 
conditions. 

Local competitive accretion occurs because the timescale for the fragmentation of subclumps into cores 
is shorter than the timescale for these clumps to merge and form the final cluster. This could possibly be 
affected by turbulence or feedback, since these processes could delay core formation inside subclumps. In 
particular, the feedback from the first core forming in a given subclump could possibly prevent the formation 
of other cores in the same subclump, so that the final cluster would form by the merger of several subclumps, 
each one containing only one core. It remains to be seen if the effect of feedback could be that drastic. 

We need a particular mass or length scale in the problem to determine the mass range of the IMF. In 
our simulations, the gas is isothermal, and there are no physical scale in the system. It is then a numerical 
scale, the minimum Jeans mass resolved by the algorithm, that determines the lower edge of the IMF. In 
simulations with a barotropic equation of state, the scale corresponds to the Jeans mass at the density where 
the gas becomes adiabatic. Turbulence and feedback could introduce a physical scale as well. But in all 
cases, there is no apparent reason for having a relationship between the mean initial Jeans mass and the 
lower edge of the IMF. Whatever the initial Jeans mass is, we expect the cloud to fragment down to the 
lowest mass scale allowed by either the physical processes or the numerical resolution. 

9. SUMMARY AND CONCLUSION 

This paper presented simulations of the fragmentation of a molecular cloud in which artificial fragmen- 
tation is prevented by using particle splitting, a technique that had not yet been applied to this kind of 
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problem. With this technique, we can follow the evolution of the entire cloud, from the initial conditions to 
the formation of a star cluster, while fully resolving the Jeans mass at all densities. 

The main objective of this paper was to demonstrate the feasibility of simulating the formation of a star 
cluster by cloud fragmentation, while properly resolving the Jeans mass throughout the entire simulation. 
We have successfully shown that, with the use of particle splitting, this can be achieved at a small fraction 
of the computational cost of a standard high-resolution or a "zoom-in" simulation. In particular, our largest 
simulation (Run C), which fully resolves the Jeans mass, has an effective resolution of 256 3 , or 17 million 
particles, while the actual number of particles in that simulation varies but never exceeds 1.1 million. The 
gain in performance is huge, and would only increase with additional splitting generations. 

We have identified four distinct phases in the evolution of the cloud, corresponding to four different 
regimes. Initially, the cloud is in the growth regime. Ovcrdense regions become denser and eventually form 
a network of filaments. Then, in the collapse regime, the gas in dense regions is converted into protostellar 
cores. In the accretion regime, the remaining gas, which started up in undcrdense regions, is mainly accreted 
by the existing cores. Finally, in the N-body regime, most of the gas has disappeared, and the evolution 
of the cluster is governed by N-body dynamics. These various regimes were certainly present in previous 
simulations, such as the ones of KB, but were not explicitly identified. The existence of a collapse and an 
accretion phase, and the fact that roughly 50% of the gas is removed during the collapse phase and 50% 
during the accretion phase, is most likely a mere consequence of the Gaussianity of the initial conditions, 
though this remains to be tested. 

In the lowest-resolution simulation, we have noticed several interesting phenomena, such as early ejec- 
tion, local competitive accretion, late starburst, and runaway accretion. Early ejections (also noticed by KB 
and Bate, Bonnell, & Bromm 2002b) occur during few-body encounters, and explain how cores that form 
early stop accreting and end up having a low mass. Local competitive accretion is a new phenomenon that 
we have identified. It occurs when several clumps of gas fragment to form sub-clusters of cores that later 
merge to form the final cluster. The first core formed in each clump does not compete for accretion until 
other cores form in the same clump, and can therefore reach high masses even if it formed late in the overall 
simulation. This explain most of the high-mass peaks at high birth ranks seen in Figure 8. Late starburst 
is caused by gas that was never dense enough to form cores until it falls into the main cluster at late times 
and gets suddenly shocked to very high densities. Runaway accretion occurs when the first core formed in 
the simulation formed significantly earlier that the next cores, and thus gets an early start in accreting gas. 
The mass difference then increases as a result of competitive accretion, and that core ends up containing a 
large fraction (up to one third) of the mass of the entire system. While all these phenomena were observed 
in the lowest-resolution simulation, only local competitive accretion was observed in the higher-resolution 
simulations. This suggests that local competitive accretion is a fundamental process that greatly affects the 
evolution of the system, while the other phenomena are less fundamental, and might not have occurred at 
all (or to a very different extent) if we had used different initial conditions. We believe that it is the small 
number of cores formed in the lowest-resolution simulation, and not its inability to resolve the Jeans mass, 
that is responsible for the occurrence of these phenomena. 

The final distribution of the core masses are roughly log-normal (in agreement with the isothermal 
simulations of KB), except for the lowest-resolution simulation, where the runaway accretion of core 
results in a very skewed distribution. The location of the distribution shifts to lower masses as the resolution 
increases, until the resolution is sufficient to resolve the Jeans mass. This is a consequence of our decision 
to follow "Scenario II" for the simulations with N gCD = 1 and 2, requiring that dense clumps are converted 
into cores only if they contain at least n m i n particles. Had we followed Scenario I, all simulations would have 
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produced (presumably) distributions with the same mean. We found that the mean of the distribution is 
determined entirely by the resolution, and not by the mean Jeans mass at the initial time, in contradiction 
with the claim of KB. 
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A. INITIAL CONDITIONS 

The technique used for generating initial conditions is a generalization of the Zel'dovich approximation 
commonly used for cosmological simulations. The same technique was used by KB, but their description 
lacks details. In this appendix, we provide a more detailed description. Notice that most of this derivation 
can be found in §10.2 of Coles & Lucchin (1995). 

To set up initial conditions, we lay down N equal-mass particles on a uniform cubic lattice inside the 
computational volume, and displace these particles to represent the density perturbation as a Gaussian 
random field with a particular power spectrum. We then adjust the particle velocities by requiring that the 
perturbation is growing with time. First, we start with equations (1), (2), and (4). We eliminate p using 
p = p(l + 5), where 5 is the density contrast. Then, assuming that the perturbation is initially small, we 
neglect the pressure gradient, and only keep terms that are linear in 8, v, and <p. Equations (1), (2), and (4) 
reduce to 

g + V.v.O, (Al) 

£— V*. (A2) 
V 2 (f> = inGpS . (A3) 

To solve these equations, we first take the divergence of equation (A2). This introduces a term in V • v, 
which we eliminate using equation (Al), and a term in V 2 0, which we eliminate using equation (A3). We 
get 

"" S =4nGp6 = Lu 2 5, (A4) 



dt 2 

where u> = {AirGp) 1 / 2 . The general solution is 

5(t, t) = A(r)e ut + B(r)e" ui . (A5) 

We assume that at the initial time U, the perturbation is in a pure growing mode, and therefore B = 0. 
Equation (A5) reduces to 

5(r,U) = A(r)e^ . (A6) 

To solve for the velocity, we decompose the density contrast and velocity into sums of plane waves, 

S(r,U) = e^]TA ke - 2 ™ k - r , (A7) 
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v(r,tO = ^v k e- 2 " k ' r , (A8) 

k 

where the sums are over all plane waves that satisfy the periodic boundary conditions and whose frequency 
does not exceed the Nyquist frequency, 

k = {n 1 ,n 2 ,n 3 )/L hox , m, n 2 , n 3 = 0, 1, 2, . . .n nyq , (A9) 

where n nyq = A 1//3 /2. We substitute equations (A7) and (A8) into equation (Al), which becomes 

U ^ke~ 27r ' k ' r = 2™ YJ k • v k e- 27rlk ' r . (A10) 



we"** 



Since the plane waves are orthogonal functions, equation (A10) implies 

= 2mk ■ v k . (All) 

To solve for v k , we make the assumption that the velocity field is vorticity-free. This implies that k x v k = 0, 
and therefore v k = |v k |k/fc and k • v k = |v k |fc, where k = |k|. Equation (All) then gives 

Vk = "^ e *' (A12) 

and equation (A8) becomes 

k 

To compute the particle displacements (Ar)j, we integrate the velocity between times — oo and ti, 

(Ar), = -— V 4#e- 2 ™ k - r = ^ • (A14) 

k 

For a given set of complex amplitudes the displacements and velocities can be computed using equa- 
tion (A13) and (A14). The amplitudes can be written as 

A k = |^k|e^ k , (A15) 

where k is a random phase with uniform probability between and 2ir, and |A k | is related to the power 
spectrum by 

P{k) = \A^\ 2 . (A16) 

Notice that the phases must satisfy the condition 0_ k = — </> k for Ar and v to be real. We assume a power 
spectrum of the form 

P(k) = ck- n , (A17) 

where c is a constant. The displacements and velocity depends upon c and ti only through the quantity 
c i/2 e ^ti_ rp n j s q uan tity i s arbitrary, but must be chosen sufficiently small that S(r, ti) <C 1, otherwise this 
linear treatment is not valid. As KB point out, it must not be chosen too small either, otherwise the early 
evolution of the system would proceed very slowly, resulting in an excessive amount of CPU time just to 
get to the nonlinear regime. We generated initial conditions by imposing that the largest particle initial 
displacements (Ar)j are equal to 10% of the mean particle spacing. This limits the initial values of S(r,ti) 
to be in the range —0.3 < 5(r, ti) < 0.3. 



-42 - 



REFERENCES 

Alecian, G., & Leorat, J. 1998, A&A, 196, 1 

Andre, P., Ward-Thompson, D., & Barsony, M. 2000, in Protostar and Planets IV, eds. V. Manning, A. 
Boss, & S. Russell (Tucson: University of Arizona Press), p. 59 

Barger, A. J., Cowie, L. L., Sanders, D. B., Fulton, E., Taniguchi, Y., Sato, Y., Kawara, K., & Okuda, H. 
1998, Nature, 394, 248 

Bate, M. R. 1998, ApJ, 508, L95 

Bate, M. R., & Bonncll, I. A. 2005, MNRAS, 356, 1201 
Bate, M. R., Bonncll, I. A., & Bromm, V. 2002a, MNRAS, 332, L65 
Bate, M. R., Bonncll, I. A., & Bromm, V. 2002b, MNRAS, 336, 705 
Bate, M. R., Bonncll, I. A., & Bromm, V. 2003, MNRAS, 339, 577 
Bate, M. R., Bonncll, I. A., & Price, N. M. 1995, MNRAS, 277, 362 
Bate, M. R., & Burkert, A. 1997, MNRAS, 288, 1060 

Bcckwith, S. V. W., Sargent, A. I., Chini, R. S., & Giisten, R. 1990, AJ, 99, 924 
Bonncll, I. A., Bate, M. R., Clarke, C. J., & Pringlc, J. E. 1997, MNRAS, 285, 201 
Bonncll, I. A., Bate, M. R., Clarke, C. J., & Pringlc, J. E. 2001, MNRAS, 323, 785 
Boss, A. P. 1998, ApJ, 501, L77 

Bromm, V., Coppi, P. S., & Larson, R. B. 1999, ApJ, 527, L5 
Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23 
Burkert, A., Bate, M. R, & Bodcnhcimcr, P. 1997, MNRAS, 289, 497 
Chapman, S. C, Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772 

Churchwcll, E. B. 1993, in Massive Stars: Their Lives in the Interstellar Medium, eds. J. P. Cassinclli & E. 
B. Churchwell (San Francisco: APS), p. 35 

Cimatti, A., Andrcani, P., Rottgering, H., & Tilanus, R. 1998, Nature, 392, 895 

Clarke, C. J., Bonncll, I. A., & Hillenbrand, L. A. 1999, in Protostar and Planets IV, eds. V. Manning, A. 
Boss, & S. Russell (Tucson: University of Arizona Press), p. 151 

Coles, P., & Lucchin, F. 1995, Cosmology. The Origin and Evolution of Cosmic Structure (New York: Wiley) 

Cox, P. et al. 2002, A&A, 387, 406 

Cumming, A., Marcy, G. W., & Butler, R. P. 1999, ApJ, 526, 890 

Dale, J. E., Bonncll, I. A., Clarke, C. J., & Bate, M. R. 2005, MNRAS, 358, 291 

Dobbs C. L., Bonncll, I. A., & Clark, P. C. 2005, MNRAS, 360, 2 



-43 - 



Doty, S. D., & Neufeld, D. A. 1997, ApJ, 489, 122 
Eisner, J. A., & Carpenter, J. M. 2003, ApJ, 598, 1341 

Elmegreen, B. G. 1985, in Protostar and Planets II, eds. V.D. C. Black, & M. S. Matthews (Tucson: 
University of Arizona Press), p. 33 

Elmegreen, B. G., Efrcmov, Y., Pudritz, R. E., & Zinnecker, H. 2000, in Protostar and Planets IV, eds. V. 
Manning, A. Boss, & S. Russell (Tucson: University of Arizona Press), p. 179 

Evans, N. J. 1999, ARA&A, 37, 311 

Evans, N. J., Rawlings, J. M. C, Shirley, Y. L., & Mundy, L. G. 2001, ApJ, 557, 193 

Evans, N. J., Shirley, Y. L., & Mueller, K. E., & Knez, C. 2002, in Hot Star Workshop III: The Earliest 
Phases of Massive Star Birth, cd. P. A. Crowther (ASP Conference Series), 267, 17 

Evrard, A. E. 1988, MNRAS, 235, 911 

Goldreich, P., & Kwan, J. 1974, ApJ, 189, 441 

Hauser, M. G. et al. 1998, ApJ, 508, 25 

Hockney, J. R., & Eastwood, J. W. 1981, Computer Simulation Using Particles (New York: McGraw-Hill) 
Jappscn, A.-K., & Klesscn, R. S. 2004, A&A, 423, 1 

Jappsen, A.-K., Klesscn, R. S., Larson, R. B., Li, Y., & Mac Low, M.-M., A&A, 435, 611 

Johnstone, D., Wilson, C. D., Moriarty-Schievcn, G., Joncas, G., Smith, G., Gregersen, E., & Fich, M. 2000, 
ApJ, 545, 327 

Jorissen, A., Mayor, M., & Udry, S. 2001, A&A, 379, 992 

Keto, E. R., Lattanzio, J. D., & Monaghan, J. J. 1991, ApJ, 383, 639 

Kitsionas, S., & Whitworth, A. P. 2002, MNRAS, 330, 129 

Klein, R. I., Fisher, R. T., McKee, C. F., & Truelovc, J. K. 1999, in Numerical Astrophysics, eds. S. M. 
Miyama, K. Tomisaka, & T. Hawana (Boston: Kluwer Academic Press), p. 131 

Klessen, R. S. 2001a, ApJ, 550, L77 

Klessen, R. S. 2001b, ApJ, 556, 837 

Klessen, R. S. 2002, private communication 

Klessen, R. S., & Burkert, A. 2000, ApJS, 128, 287 

Klessen, R. S., & Burkert, A. 2001, ApJ, 549, 386 

Klessen, R. S., Burkert, A., & Bate, M.R. 1998, ApJ, 501, L205 

Klessen, R. S., Heitsch, F., & Mac Low, M.-M. 2000, ApJ, 535, 887 

Kurtz, C, Cesaroni, R., Churchwell, E. B., Hofncr, P., & Walmsley, M. 2000, in Protostar and Planets IV, 
eds. V. Manning, A. Boss, & S. Russell (Tucson: University of Arizona Press), p. 299 



-44- 



Lada, E. A., Evans, N. J., & Falgarone, E. 1997, ApJ, 488, 286 
Lada, E. A., & Lada, C. J. 1995, AJ, 109, 1682 
Lada, C. J. & Lada, E. A. 2003, ARA&A, 41, 57 

Lanzetta, K. M., Yahata, N., Pascarelle, S., Chen, H., & Fernandez-Soto, A. 2002, ApJ, 570, 492 
Larson, R. B. 1978, MNRAS, 184, 69 

Larson, R. B. 2003, in Galactic Star Formation Across the Stallar Mass Spectrum, eds. J. M. De Vuizer & 
N. S. van der Bliek, ASP Conference Series Vol. 287, p. 65. 

Leung, C. M. 1975, ApJ, 199, 340 

Li, Y., Klcsscn, R. S., & Mac Low, M.-M. 2003, ApJ, 592, 975 

Mac Low, M. & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125 

Madau, P. 1999, in After the Dark Ages: When Galaxies were Young (the Universe at 2 < z < 5), eds. S. 
Holt & E. Smith (AIP Press), p. 299 

Marcy, G., Butler, R. P., Fischer, D., Vogt, S., Wright, J. T., Tinncy, C. C, & Jones, H. R. A. 2005, Progress 
of Theoretical Physics Supplement, 158, 24 

McKee, C. F., & Tan, J. C. 2002, Nature, 416, 59 

McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850 

McWilliam, A. 1997, ARA&A, 35, 503 

Monaghan, J. J. 1992, ARA&A, 30, 543 

Mottc, F., Andre, P., & Ncri, R. 1998, A&A, 336, 150 

Mueller, K. E., Shirley, Y. L., Evans, N. J., & Jacobson, H. R. 2002, ApJS, 143, 469 

Mundy, L. G., Looney, L. W., & Welch, W. J. 2000, in Protostar and Planets IV, eds. V. Manning, A. Boss, 
& S. Russell (Tucson: University of Arizona Press), p. 355 

Myers, P. C, Evans, N. J., & Ohashi, N. 2000, in Protostar and Planets IV, eds. V. Manning, A. Boss, & S. 
Russell (Tucson: University of Arizona Press), p. 217 

O'Dell, C. R., & Wen, Z. 1994, ApJ, 436, 194 

Ossenkopf, V. & Hcnning, T. 1994, A&A, 291, 943 

Owen, J. M., Villumscn, J. V., Shapiro, P. R., & Martcl, H. 1998, ApJS, 116, 155 
Padoan, P., & Nordlund, A. 2002, ApJ, 576, 870 

Plume, R., Jaffe, D. T., Evans, N. J., Martin-Pintado, J., & Gomez-Gonzalez, J. 1997, ApJ, 476, 730 

Pudritz, R. E., & Fiegc, J. D. 1999, in New Perspectives on the Interstellar Medium, eds. A. R. Taylor, T. 
L. Landecker, & G. Joncas (San Francisco: APS), p. 235 



-45 - 



Puget, J.-L., Abergal, A., Bernard, J. -P., Boulanger, F., Burton, W. B., Desert, F. X., & Hartmann, D. 
1996, A&A, 308, L5 

Scalo, J., Vazquez-Semadeni, E., Chappcll, D., & Passot, T. 1998, ApJ, 504, 835 
Schmeja, S., & Klesscn, R. S. 2004, A&A, 419, 405 

Shapiro, P. R., Martel, H., Villumscn, J. V., & Owen, J. M. 1996, ApJS, 103, 269 

Shirley, Y. L., Evans, N. J., Young, K. E., Knez, C, & Jaffe, D. T. 2003a, ApJS, 149, 375 

Shirley, Y. L., Mueller, K. E., Young, C. H., & Evans, N. J. 2003b, ASP Conf. Ser. 287: Galactic Star 
Formation Across the Stellar Mass Spectrum, 298 

Schneider, S., & Elmcgrccn, B. G. 1979, ApJS, 41, 87 

Shu, F. H. 1977, ApJ, 214, 488 

Shu, F. H., Adams, F. C, & Lizano, S. 1987, ARA&A, 25, 23 

Smail, I., Ivison, R., Blain, A., Kneib, J. -P., & Owen, F. 2000, in Imaging the Universe in Three Dimensions: 
Astrophysics with Advanced Mulit- Wavelength Imaging Devices, eds. W. van Breugel & J. Bland- 
Hawthorne (San- Francisco: ASP), p. 248 

Strom, S. 1995, Rev.MexA.&A., 1, 317 

Terebey, S., Shu, F. H., & Casscn, P. 1984, ApJ, 286, 529 

Testi, L., & Sargent, A. I. 1998, ApJ, 508, L91 

Tillcy, D. A., & Pudritz, R. E. 2004, MNRAS, 353, 769 

Tohline, J. E. 1982, Fundam. Cosmic Phys., 8, 1 

Truelove, J. K., Klein, R. I., McKee, C. F., Holliman, J. H., Howell, L. H., & Greenough, J. A. 1997, ApJ, 
489, L179 

Walter, F., Carilli, C, Bcrtoldi, F., Mcnten, K., Cox, P., Lo, K. Y., Fan, X., & Strauss, M. A. 2004, ApJ, 
615, L17 

Wang, Y., Jaffe, D. T., Evans, N. J., Hayashi, M., Tatcmatsu, K., & Zhou, S. 1993, ApJ, 419, 707 

Young, K. E., Lee, J., Evans, N. J., Goldsmith, P. F., & Doty, S. D. 2004, ApJ, 614, 252 

Zhou, S., & Evans, N. J. 1994, in Clouds, Cores, and Low-Mass Stars, eds. D. P. Clemens & R. Barvainis 
(San Francisco: APS), p. 183 

Zhou, S., Evans, N. J., Kompe, C, & Walmsley, C. M. 1993, ApJ, 404, 232 



This preprint was prepared with the AAS IATjTjX macros v5.2. 



